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1.  Introduction 

Computational  Fluid  Dynamics  (CFD)  has  enormous  potential  to  impact  the  analysis,  design,  and 
optimization  of  U.S.  Air  Force  systems.  However,  the  predictive  capability  of  CFD  depends  not  only  on 
the  validity  of  the  sub-models  employed  (e.g.,  turbulence,  chemistry,  multi-phase  flow)  and  the 
uncertainties  present  in  the  system  and  surroundings,  but  also  on  the  ability  to  reliably  estimate  and 
reduce  numerical  errors  (Oberkampf  and  Roy,  2010;  Roy  and  Oberkampf,  201 1).  While  there  are  many 
sources  of  numerical  error  in  a  CFD  simulation,  the  largest  and  most  difficult  source  to  estimate  (and  that 
our  proposed  effort  targets)  is  the  error  related  to  the  resolution  and  quality  of  the  spatial  mesh,  i.e.,  the 
spatial  discretization  error.  For  example,  the  recent  AIAA  Drag  Prediction  Workshops  examined 
simplified  transport  aircraft  challenge  cases  computed  by  numerous  CFD  practitioners;  however,  even 
after  three  such  workshops,  analysis  of  the  results  found  that  “grids  remain  a  first  order  effect”  (Morrison 
and  Hemsch,  2007).  While  this  proposal  is  focused  on  the  reliable  estimation  and  control  of  discretization 
error  in  CFD,  the  proposed  techniques  also  apply  to  other  areas  such  as  computational  structural 
mechanics,  computational  heat  transfer,  etc. 

Most  compressible  CFD  codes  employ  finite  volume  or  finite  difference  discretizations.  While  there 
has  been  a  great  deal  of  work  on  discretization  error  estimation  for  finite  element  method  (mainly  for 
elliptic  problems),  compressible  CFD  codes  tend  to  rely  on  simple  approaches  such  as  Richardson 
extrapolation.  The  main  drawback  to  Richardson  extrapolation  is  that  it  requires  two,  or  even  three, 
systematically-refined  grids  in  order  to  obtain  the  error  estimate.  For  most  practical  applications,  each 
grid  cell  is  usually  refined  by  a  factor  of  two  for  each  spatial  direction,  resulting  in  eight  new  fine  grid 
cells  being  generated  for  each  coarse  grid  cell.  For  3D  applications,  if  the  coarse  grid  uses  10  million 
cells,  then  two  additional  levels  of  refinement  would  involve  80  million  and  640  million  cells.  Such  large 
cell  counts  are  prohibitively  expensive,  especially  when  the  original  coarse  grid  provides  adequate 
solutions  and  the  additional  grid  levels  are  simply  required  for  the  error  estimates.  Reliable  discretization 
error  estimates  are  needed  both  for  the  estimation  of  total  predictive  uncertainty  in  a  simulation  and  to 
provide  a  stopping  criteria  for  mesh  adaptation  strategies. 

Automatic  mesh  adaptation  for  CFD  has  been  a  goal  of  researchers  for  more  than  three  decades. 
However,  an  examination  of  current  government  and  industry  CFD  codes  shows  that  the  dream  of  robust 
and  automatic  mesh  adaptation  has  not  yet  been  realized.  The  few  CFD  codes  that  do  claim  to  do  mesh 
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adaptation  usually  drive  it  with  simple  eriteria  sueh  as  solution  gradients  or  eurvature  (i.e.,  Hessians); 
however,  these  approaehes  ean  fail  disastrously  for  eomplex  problems  (Ainsworth  and  Oden,  2000; 
Dwight,  2008).  Mathematieally  rigorous  approaehes  for  driving  mesh  adaptation  are  diseussed  in  this 
proposal.  Onee  the  strategy  for  driving  the  adaption  is  ehosen,  a  teehnique  must  be  ehosen  for  aetually 
aeeomplishing  the  adaptation.  Adaptation  approaehes  for  eompressible  CFD  eomputations  are  praetieally 
limited  to  mesh  movement  (i.e.,  r-adaptation)  for  struetured  grids  and  r-  and  h-  adaptation  (i.e.,  loeal  eell 
additions/deletions)  for  unstruetured  grids.  Adaptation  based  on  ehanging  the  loeal  formal  order  of 
aeeuraey  of  the  seheme  (i.e.,  p-adaptation)  for  eompressible  CFD  eodes  (espeeially  on  unstruetured  grids) 
is  diffieult,  although  some  progress  is  being  made  through  the  use  of  diseontinuous  Galerkin  methods 
(Hesthaven  and  Warburton,  2008).  Automatie  mesh  adaptation  is  required  for  eomplex  CFD  problems  in 
order  to  obtain  reliable  diseretization  error  estimates. 

Inaeeurate  predietion  of  aerodynamie  loads  was  found  to  be  the  leading  eause  of  unantieipated 
struetural  response  or  damage  (Love  et  al.,  2003),  eosting  the  U.S.  Air  Foree  billions  of  dollars  and 
adversely  impaeting  system  readiness.  Our  approaeh  will  provide  designers  and  analysts  with  teehniques 
to  quantify  and  reduee  diseretization  error  in  CFD  predietions  without  the  overhead  of  ereating,  and 
solving  on,  multiple  meshes.  Our  work  will  lead  to  signifieant  improvements  in  the  aeeuraey  and 
effieieney  of  CFD  predietions  of  aerodynamie  loads.  The  outeomes  of  the  proposed  work  will  be 
improved  aerodynamie  predietions  for  preliminary  design,  parametrie  studies,  sensitivity  analysis, 
uneertainty  quantifieation,  optimization,  optimization  under  uneertainty,  and  multiphysies  eomputations. 


II.  Background 

The  diseretization  error  is  the  numerieal  error  due  to  the  mesh  and  time  step  ehosen  for  the 
simulation.  It  is  formally  defined  as  the  differenee  between  the  exaet  solution  to  the  diserete  equations 
and  the  exaet  solution  to  the  underlying  partial  differential  or  integral  equations  (referred  to  eolleetively 
as  the  PDFs  in  this  proposal).  This  relationship  is  given  in  equation  form  as 

~  (1) 


where  is  the  loeal  diseretization  error,  Uh  is  the  numerieal  solution  (negleeting  iterative  and  round-off 
error),  and  u  denotes  the  exaet  solution  to  the  PDFs. 


Truncation  Error  Framework 

The  truneation  error  is  the  differenee  between  the  diserete  equation  and  the  underlying  PDF  and  ean 
be  found  by  using  Taylor  series  expansions.  A  simple  example  of  a  truneation  error  analysis  follows. 
Consider  ID  steady  Burgers’  equation 


du  d^u 

- V - y 

dx  dx 


=  0 


(2) 


whieh  provides  a  simple,  sealar,  nonlinear  model  equation  for  the  Navier-Stokes  equations  as  it  eontains 
nonlinear  eonveetion  and  linear  diffusion  with  eonstant  viseosity  v.  Burgers’  equation  will  be  used  as  a 
simple  example  throughout  this  proposal.  For  a  simple  seeond-order  aeeurate  finite  differenee 
diseretization  of  Burgers’  equation  using  a  eonstant  mesh  spaeing  h  =  Ax,  a  truneation  error  analysis  gives 
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+  u, 


d^u 


'  dx^ 


Ax  d  u 

6  dx 


.I2 


■  + 


o[aj'] 


Th(u) 


where  Lh(u)  is  the  diserete  operator,  L(u)  the  PDE  operator  (here  it  is  restrieted  to  the  nodes),  and  Th(u)  is 
the  truneation  error.  The  diserete  equations  are  exaetly  satisfied  by  the  diserete  solution  at  the  nodes  Uh 
(i.e.,  =  0),  the  PDE  is  exaetly  solved  by  its  eontinuous  solution  (i.e.,  L{u)  =  0 ),  and  the  truneation 

error  operates  on  a  eontinuous  funetion  but  produees  diserete  values.  This  is  a  general  relationship  for 
arbitrary  (but  smooth)  u(x)  whieh  assumes  appropriate  operators  are  used  to  restriet  eontinuous  funetions 
to  diserete  points  or  prolong  diserete  funetions  to  a  eontinuous  spaee  (see  the  Teehnieal  Approaeh  seetion 
for  a  more  rigorous  diseussion).  Using  the  above  operator  notation,  this  relationship  ean  be  expressed  as 
the  Generalized  Truncation  Error  Expression  (GTEE)  (Roy,  2009  and  2010a): 


C{u)  =  L(u)  +  t^{u).  (3) 

The  relationship  between  the  truneation  error  and  the  diseretization  error  ean  be  established  by 
examining  the  eontinuous  diseretization  error  transport  equation.  Inserting  the  exaet  solution  to  the 
diserete  equations  m  into  Equation  (3),  then  subtraeting  the  original  (eontinuous)  governing  equation 
L{u)  =  0  gives 

L(Mj-L(M)  +  r^(Mj  =  0. 

If  the  operator  is  linear  (or  has  been  linearized  by  assuming  u  then  L(w^ )  -  L(u)  =  L{u^  -  u) . 

Employing  the  definition  of  the  diseretization  error  from  Equation  (1)  results  in 


(4) 


whieh  for  Burgers’  equation  beeomes 


ds.  d^s.  d^u.  d^u.  tsx^  ^ 

U  —  -V - :^  =  -U - f - +  V - f - +  0 


dx  dx 


dr  6 


dx"  12 


Ax' 


(5) 


Equations  (4)  and  (5)  represent  a  eontinuous  transport  equation  (similar  to  the  error  equation  in  finite 
elements)  whieh  ean  be  solved  on  the  same  grid  as  the  original  problem  to  estimate  the  diseretization 
error.  It  shows  that  diseretization  error  is  transported  in  the  same  manner  as  the  solution,  and  that  it  is 
loeally  generated  by  the  truneation  error.  When  the  truneation  error  is  redueed,  then  the  loeal  souree  for 
the  diseretization  error  is  also  redueed,  leading  to  less  diseretization  error  ‘‘produetion”  over  the  domain. 
Broadly  speaking,  when  globally  ‘‘good”  numerieal  solutions  are  desired,  then  the  truneation  error  is  the 
ideal  driver  for  mesh  adaptation  (Baker,  1997;  Roy,  2009).  When  only  a  small  number  of  solution 
funetionals  are  of  interest  (e.g.,  lift  or  drag  in  an  aerodynamies  simulation),  then  adjoint  methods  ean  be 
used  to  aeeount  for  the  sensitivity  of  the  solution  funetional  to  loeal  truneation  errors  (e.g.,  Pieree  and 
Giles,  2000;  Venditti  and  Darmofal,  2003). 
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Discretization  Error  Estimators 

The  discretization  error  Sh  was  given  in  Equation  (1)  and  is  formally  defined  as  the  difference 
between  the  exact  solution  to  the  discrete  equations  m  and  the  exact  solution  to  the  underlying  PDEs  u  . 
The  spatial  discretization  error  is  generally  the  largest  contributor  to  the  overall  numerical  error,  with  the 
other  contributors  being  temporal  discretization  error,  iterative  convergence  error,  and  round-off  error. 
There  are  two  types  of  discretization  error  estimators  (Roy,  2010a).  The  first  involves  comparing  the 
numerical  solution  (or  its  gradients)  to  higher-order  accurate  estimates  of  the  solution  (gradients). 
Included  in  this  category  are  Richardson  extrapolation  for  the  solution  and  solution  functionals  which 
requires  solutions  on  multiple  meshes  (Roy,  2005)  and  recovery-based  methods  for  gradients  from  finite 
elements  (e.g.,  Zienkiewicz  and  Zhu,  1992).  The  second  class  of  methods  for  estimating  discretization 
error  employs  either  the  continuous  or  discrete  residual  (not  to  be  confused  with  the  iterative  residuals 
which  are  often  used  to  monitor  iterative  convergence)  or,  equivalently,  the  truncation  error.  Residual- 
based  methods  employ  information  about  the  problem  being  solved  and  include  error  transport  equations 
(e.g.,  Cavallo  and  Sinha,  2007),  defect  correction  (e.g.,  Skeel,  1986),  adjoint  methods  (e.g.,  Venditti  and 
Darmofal,  2003),  and  explicit/implicit  residual  methods  from  finite  elements  (Ainsworth  and  Oden, 
2000).  Finally,  there  is  a  hybrid  approach  called  least  squares  extrapolation  where  the  extrapolation 
employs  unknown,  spatially  varying  coefficients  which  are  determined  by  minimizing  the  residual  on  a 
finer  mesh  (Garbey  and  Shyy,  2003). 

Spatial  discretization  error  is  the  most  difficult  type  of  numerical  error  to  estimate  reliably. 

Regardless  of  the  approach  used  to  estimate  the  discretization  error,  the  numerical  solution(s)  must  be  in 
the  asymptotic  range  in  order  to  obtain  reliable  estimates.  The  asymptotic  range  is  defined  as  the  range  of 
mesh  resolutions  where  the  discretization  error  reduces  at  the  theoretical  (i.e.,  formal)  rate  with  mesh 
refinement.  This  asymptotic  range  is  defined  only  in  terms  of  systematic  mesh  refinement  (Oberkampf 
and  Roy,  2010)  where  the  mesh  is  refined  uniformly  by  a  factor  h  in  each  coordinate  direction,  e.g.. 


Ax  Ay  Az 

^ref  ^Jref  ^ref 


(6) 


and  the  mesh  quality  is  constant  or  improves  with  mesh  refinement.  Ensuring  systematic  mesh  refinement 
can  be  challenging,  especially  for  unstructured  meshes  which  contain  more  than  one  type  of  mesh 
topology  (e.g.,  hexahedral,  tetrahedral,  and  prismatic  elements). 

Verifying  that  the  solutions  are  in  the  asymptotic  range  is  traditionally  done  by  computing  the 
observed  order  of  accuracy  using  numerical  solutions  on  three  systematically -refined  meshes.  For 
systematic  refinement  by  the  factor  r,  one  has  h/me  =  h,  hmedium  =  rh,  and  hcoarse  =  r^h  and  the  observed 
order  of  accuracy  p  can  be  found  as  (Roache,  2009): 


In 


1,2,  -I 


rh 


frh  -  fl 


\  J  rh 


h  J 


ln(r) 


(7) 


The  grid  refinement  factor  can  generally  be  as  small  r=\.\  without  round-off  and  iterative  error  polluting 
the  results;  however,  from  a  practical  standpoint,  a  refinement  factor  of  r  =  2  is  often  used  which  implies 
insertion  of  new  nodes  between  each  existing  node  for  structured  grids  and  local  sub-division  of  elements 
in  unstructured  grids.  For  three-dimensional  flows,  r  =  2  results  in  a  factor  of  8  increase  in  total  grid  cells 
for  each  level  of  refinement  and  can  thus  be  prohibitively  expensive.  Furthermore,  the  observed  order  of 
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accuracy  will  only  match  the  formal  order  when  all  three  grid  levels  are  in  the  asymptotic  range.  A  major 
challenge  is  to  obtain  reliable  discretization  error  estimates  without  requiring  solutions  with  hundreds  of 
millions  of  cells.  For  practical  problems,  it  is  our  opinion  that  the  asymptotic  range  is  not  achievable  (or  at 
least  demonstrable)  without  effective  local  mesh  adaptation  strategies. 

Truncation  Error  and  Residuals 

The  truncation  error  can  be  related  to  the  residuals  used  in  residual -based  error  estimators  and  adjoint 
methods.  Inserting  the  numerical  solution  m  into  the  GTEE  given  in  Equation  (3)  gives  the  continuous 
residual 


Liu,)  =  -T,iuJ  (8) 

which  is  analogous  to  the  finite  element  residual  (Ainsworth  and  Oden,  2000)  and  the  source  term  used  in 
the  differential  form  of  defect  correction  (Skeel,  1986).  If  instead  the  exact  solution  to  the  continuous 
governing  equation  u  (or  an  approximation  thereof)  is  inserted  into  the  GTEE,  one  obtains  the  discrete 
residual 


Lh(u)  =  Tf,iu).  (9) 

It  is  this  discrete  residual  (or  its  approximation)  that  is  approximated  in  many  adjoint  methods  (Venditti 
and  Darmofal,  2003)  and  discretization  error  transport  equations  (e.g.,  Shih  and  Williams,  2009). 

Mesh  Adaptation 

There  are  two  main  approaches  for  performing  mesh  adaption.  Mesh  refinement  (h-adaptation) 
provides  for  the  sub-division  of  mesh  cells  to  improve  mesh  resolution,  while  mesh  movement  (r- 
adaptation)  seeks  to  move  mesh  cells  from  one  region  to  another.  For  general  unstructured  grid  methods 
the  h-adaption  approach  is  the  most  popular,  while  for  structured  grid  methods  the  r-adaption  approach  is 
most  often  used.  In  both  h-  and  r-adaptation,  weighting  functions  are  generally  employed  to  drive  the 
adaptation  process.  Regardless  of  how  the  adaptation  is  performed,  the  more  difficult  challenge  is  finding 
an  appropriate  criterion  to  adapt  on.  The  mesh  adaptation  criterion  that  is  found  in  most  commercial  CFD 
codes  is  based  on  solution  features  such  as  solution  gradients,  solution  curvature,  vortex  cores,  or  shock 
waves.  However,  feature -based  adaptation  is  certainly  not  optimal  (Roy,  2009)  and  can  even  fail  to  reduce 
the  discretization  error  in  some  cases  (e.g.,  Ainsworth  and  Oden,  2000;  Dwight,  2008). 

Since  it  is  the  discretization  error  that  one  would  like  to  reduce,  at  first  glance,  it  might  appear  that  the 
discretization  error  would  be  a  good  criterion  to  drive  the  adaptation  process.  However  adaptation  based 
on  the  discretization  error  is  prone  to  adapting  to  components  of  the  discretization  error  that  have  been 
transported  from  other  regions  of  the  domain  rather  than  the  error  that  is  locally  generated  (Gu  and  Shih, 
2001;  Roy,  2009).  The  successful  use  of  recovery-based  mesh  adaptation  in  finite  elements,  where 
adaptation  is  driven  by  the  estimated  discretization  error  in  the  solution  gradients,  relies  on  the 
superconvergence  property  of  finite  element  methods  (Ainsworth  and  Oden,  2000),  as  well  as  possibly  the 
elliptic  mathematical  character  of  the  equations.  A  recovery-type  method  for  finite  volume  schemes  has 
also  been  developed  based  on  the  solution  interpolation  error  (e.g.,  see  McRae  2000).  Recovery-based 
methods  are  not  appropriate  for  driving  mesh  adaptation  in  cases  where  the  governing  equations  are 
hyperbolic  or  when  finite  difference  or  finite  volume  methods  are  used. 
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Examination  of  the  continuous  discretization  error  transport  equation  (called  the  error  equation  in 
finite  elements)  given  by  Equation  (4)  shows  that  underlying  PDE  transports  the  discretization  error  in  the 
same  fashion  as  the  solution  properties.  For  our  Burgers’  equation  example,  the  error  transport  equation 
given  by  Equation  (5)  shows  that  the  error  will  be  both  convected  at  the  local  velocity  and  locally  diffused 
according  to  the  viscosity.  More  importantly,  it  shows  that  the  truncation  error  serves  as  a  local  source 
term  for  the  discretization  error;  thus  it  is  the  truncation  error  (or  its  approximation)  that  should  be  used  to 
drive  the  mesh  adaptation  process.  The  simplest  approach  is  to  use  the  magnitude  of  the  local  truncation 
error  to  drive  the  adaptation  process.  A  more  advanced  approach  is  to  account  for  both  the  mesh 
resolution  components  and  the  mesh  quality  components  (stretching,  skewness,  aspect  ratio,  etc.)  when 
adapting  (Yamaleev,  2001;  Alyanak  et  ah,  201 1).  Finally,  if  one  is  interested  in  only  a  small  number  of 
solution  functionals  (e.g.,  lift  and  drag),  adjoint  methods  can  be  used  to  drive  adaption  by  including 
sensitivities  of  those  functionals  to  the  local  truncation  error  (e.g.,  Venditti  and  Darmofal,  2003;  Wang 
and  Mavriplis,  2009;  Fidkowski  and  Darmofal,  201 1). 

III.  Objective 

The  overall  objective  of  our  effort  was  to  develop  and  demonstrate  robust,  reliable,  and  automatic 
methods  for  controlling  discretization  error  for  CFD  applications.  These  methods  are  based  on  rigorous 
theory  rather  than  ad  hoc  approaches  such  as  feature -based  adaptation.  Given  the  ambitious  nature  of  this 
objective,  we  limit  the  scope  of  the  problem  by  restricting  our  applications  to  steady-state  CFD  problems 
on  structured  grids;  however,  we  specifically  targeted  approaches  which  are  extendible  to  unstructured 
grids.  The  focus  is  on  practical  CFD  applications,  with  the  ultimate  application  being  the  steady  3D 
Reynolds-Averaged  Navier-Stokes  (RANS)  equations.  The  RANS  equations  are  chosen  for  their 
applicability  to  design,  analysis,  and  optimization  of  aerospace  flight  and  propulsion  systems.  The 
proposed  approaches  are  applicable  to  all  discretization  methods  (finite  difference,  finite  volume,  and 
finite  element)  and  ultimately  extendable  across  the  entire  Mach  number  range  including  subsonic, 
transonic,  supersonic,  and  hypersonic  speeds. 

IV.  CFD  Code 

The  CFD  code  that  served  as  the  basis  for  this  work  is  the  finite  volume  Euler/Navier-Stokes  solver 
developed  by  our  group  (Derlaga  et  al.,  2013).  The  code  employs  structured  grids  and  has  options  for 
various  explicit  and  implicit  iterative  schemes.  Several  code  enhancements  were  made  during  this  effort 
including  extension  to  laminar  Navier-Stokes,  3D,  and  parallelism.  All  coding  implementations  were 
done  in  a  modular  fashion  to  allow  future  code  development. 

V.  Technical  Approach 

Our  work  can  be  grouped  into  four  broad  tasks.  Our  initial  efforts  focused  on  methods  for  estimating 
the  truncation  error.  The  second  task  addressed  the  further  development  and  refinement  of  the  GTEE 
framework  including  the  required  properties  of  the  interpolation  operators  to  be  used.  The  third  task 
involved  the  development  and  evaluation  of  residual-based  discretization  error  estimators  within  the 
GTEE  framework  and  includes  error  transport  equation,  defect  correction,  and  adjoint  methods.  The  final 
task  was  to  develop  and  evaluate  mesh  adaptation  strategies  that  are  based  on  the  residual/truncation 
error. 
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VI.  Truncation  Error  Estimation 


Framework  for  Residual-Based  Discretization  Error  Estimation 

As  discussed  earlier,  the  original  Generalized  Truncation  Error  Expression  (GTEE)  requires 
appropriate  operators  to  prolong  discrete  quantities  (e.g.,  finite  volume  solutions  at  cell  centers)  to 
continuous  functions  and  to  restrict  continuous  functions  to  cell  center  or  nodal  locations  (Roy,  2009).  We 
introduce  the  interpolation  function  1  which  can  perform  both  prolongation  and  restriction  operations. 
This  interpolation  function  is  designed  to  be  read  from  bottom  to  top.  Consider  the  following  examples 
where  is  a  discrete  numerical  solution  on  a  fine  mesh,  ^2/^  is  a  discrete  numerical  solution  on  a  coarse 

mesh,  and  u  is  the  continuous  exact  solution  to  the  governing  PDEs: 

•  prolongation  of  to  a  continuous  space:  I 

•  prolongation  of  to  mesh  h:  ^ih^ih 

•  restriction  of  to  mesh  2h: 

•  restriction  of  u  to  mesh  h\  I^u 

When  no  subscript  or  superscript  is  present,  a  continuous  function  is  assumed.  Using  this  interpolation 
operator,  the  GTEE  can  be  rigorously  recast  for  finite  difference  and  finite  volume  schemes  as 

4(/'-w)  =  /^L(w)  +  r,(w)  (10) 

for  any  general  smooth  function  w  .  The  three  terms  in  Equation  (10)  refer  to  nodal  values  for  finite 
difference  methods  and  cell-averaged  values  for  finite  volume  methods.  In  addition,  the  discretization 
error  definition  from  Equation  (1)  can  be  formulated  either  continuously  or  discretely  as 

£  =  I^^U^-U  or  £^=U^-I'‘u,  (11) 


respectively. 

If  we  consider  a  discretized  function,  then  the  GTEE  can  be  written  as 

4K)  =  ^'L{lh^h)  +  hihr^h)  ■  (12) 

The  prolongation  from  discrete  space  to  continuous  space  is  an  approximation  that  is  better  written  as 

/,=/;+o(/.«*')  (13) 

to  account  for  the  error  in  the  prolongation,  where  q  is  the  order  of  the  reconstructed  polynomial  for 
which  the  k-exact  and  ENO  schemes  are  q+l  order  accurate.  Equation  (12)  is  then  written  as 

4(w.)  =  /*L(/>.+0(F*'))+r.(/;«..)+0(A>*').  ,14) 


or 
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(15) 


L,xwj=i''L(iiw,)+T,(nw,)+om. 

where  q  is  an  order  of  aeeuraey  whieh  is  related  to  q  but  may  be  modified  by  nonlinearities  in  the  L(  j 

and  Th(  j  operators.  Due  to  the  error  in  the  prolongation,  it  eannot  be  generally  assumed  that  Wh  =  I^hwh 
unless  finite  volume  eonsistent  reeonstruetions/prolongations  are  used,  i.e.,  reeonstruetions  that  will 
exaetly  integrate  out  to  the  average  value  over  the  eell. 

As  is  evident  from  the  GTEE,  the  truneation  error  is  the  differenee  between  the  diserete  and 
eontinuous  governing  equations  and  has  been  shown  to  be  the  loeal  souree  of  diseretization  error  in  a 
numerieal  solution  (Roy,  2009)  where  the  diserete  form  of  the  diseretization  error  is  defined  as  the  right 
side  of  Equation  (1 1).  Here,  m  is  the  exaet  solution  to  the  diserete  equations  sueh  that  Lh(uh)  =  0,  and  u  is 
the  exaet  solution  to  the  original  governing  equations  sueh  that  L(u)  =  0  .  The  aeeurate  estimation  of 
truneation  error  is  a  key  aspeet  to  several  residual -based  diseretization  error  estimation  methods  sueh  as 
the  error  transport  equations  (Zhang  et  ah,  2002;  Qin  and  Shih,  2002),  defeet  eorreetion  (Pereyra,  1965; 
Stetter,  1978;  Skeel,  1986),  and  adjoint  methods  (Giles  and  Pieree,  2000).  Also,  as  the  loeal  souree  of 
diseretization  error,  truneation  error  has  been  shown  to  be  an  effeetive  mesh  adaption  indieator  (Roy, 
2009),  possibly  with  adjoint  weighting  (Venditti  and  Darmofal,  2003). 

The  foeus  of  this  report  is  a  truneation  error  estimation  method  whieh  requires  a  prolongation  of  the 
numerieal  solution  to  a  eontinuous  spaee.  Inserting  the  exaet  solution  to  the  diserete  equations  m  into 
Equation  (15),  and  noting  that  Lh(uh)  =  0,  results  in 

h(i>,)=-i'-L(nu„)+om.  (IS, 

The  aeeuraey  of  the  truneation  error  estimate  is  determined  by  eomparing  the  estimated  truneation  error  to 
the  exaet  truneation  error  whieh  requires  the  exaet  solution  to  the  governing  equations.  It  ean  be  found  by 
inserting  the  exaet  solution  into  the  GTEE  and,  noting  that  L(u)  =  0  ,  gives 

where  the  exaet  and  approximate  truneation  errors  are  related  by 

Various  method  of  estimating  truneation  error  have  been  developed  and  tend  to  be  speeifie  to  the 
applieation.  For  example,  adjoint  methods  typieally  use  an  embedded  grid  approaeh  developed  by 
Venditti  and  Darmofal  (2003).  This  method  inserts  a  eoarse  grid  solution  into  a  finer  grid  using  a 
reeonstruetion  method 

h(W=-iLh,AC’M+o(h'). 

where  h  represents  the  eoarse  grid  and  h/r  represents  the  fine  grid  (refined  by  faetor  r).  We  modified  this 
method  (Phillips  and  Roy,  201 1)  to  improve  aeeuraey  sinee  the  truneation  error  on  the  fine  grid  solution 
is  still  signifieant  for  seeond-order  sehemes  with  typieal  refinement  faetors  of  two.  Equation  (19)  is 


(17) 

(18) 

(19) 
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modified  to  include  the  effects  the  fine  grid  truncation  error  assuming  that  the  truncation  error  is 

h 


asymptotic  so  that  where  pf  is  the  formal  order  of  accuracy  of  the  truncation  error: 


rPf 


.Pf 


-1 


+  0{h^). 


(20) 


A  detailed  derivation  is  given  in  by  Phillips  and  Roy  (201 1).  The  formal  order  of  truneation  error  is 
related  to  the  formal  order  of  aeeuraey  of  the  diseretization  error;  however,  for  unstruetured  grids  or  grids 
with  non-smooth  transformations  (e.g.,  with  randomly  perturbed  nodes),  the  formal  order  of  truneation 
error  ean  be  lower  than  the  formal  order  of  the  diseretization  error  (Diskin  and  Thomas,  2007).  Shih  and 
Williams  (2009)  estimated  truneation  error  for  the  error  transport  equations  by  inserting  a  finer  grid 
solution  into  the  eoarse  grid  diserete  operator 


~  ^hih/r^h/r)  * 


(21) 


Fulton  (2003)  used  a  similar  embedded  grid  approaeh  for  finite  differenee  diseretization  sehemes,  but 
adjusted  the  truneation  error  estimate  to  take  into  aeeount  the  relative  magnitudes  of  the  truneation  error 
between  the  two  grids 


~  ^hihlr^hlr) 

V 


Pf 


(22) 


Fraysse  et  al.  (2012)  extended  Fulton’s  work  to  finite  volume  methods  and  ineluded  the  effeets  of 
iterative  eonvergenee  error  for  numerieal  solutions  whieh  are  not  fully  eonverged.  Under  this  researeh 
grant,  we  followed  a  similar  approaeh  but  re-interpolated  the  estimated  truneation  error  baek  onto  the 
eomputational  grid 


1 


1 


(23) 


A  detailed  derivation  is  given  in  by  Phillips  and  Roy  (201 1). 


Solution  Reconstruction 

The  eell  average  of  the  reeonstrueted  solution  is  eomputed  using  a  Curtis-Clenshaw  (Clenshaw  and 
Curtis,  1960)  quadrature  whieh  is  extended  to  multiple  dimensions  using  Smolyak’s  (Smolyak,  1963) 
sparse  grid  eonstruetion.  We  ehose  Curtis-Clenshaw  over  the  more  eommon  Gauss  quadrature  beeause 
the  integration  points  are  nested  (i.e.,  higher-order  reeonstruetions  use  all  of  the  lower  order  quadrature 
points)  thus  allowing  adaptive  reeonstruetion  and  inexpensive  quadrature  error  estimates.  The  Curtis- 
Clenshaw  quadrature  exaetly  integrates  a  ^-th  order  polynomial  using  q+l  funetion  evaluations  in  one 

dimension.  The  quadrature  domain  is  a  hypereube  over  the  range  0<  ^  <  1  so  for  integration  of 
arbitrary  eells  a  Ist-order  polynomial  is  used  to  map  the  arbitrary  eell  to  the  quadrature  domain  where  the 
highest  polynomial  order  for  eaeh  dimension  is  one.  For  example,  in  two  dimensions  the  polynomial  takes 
the  form 
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(24) 


This  polynomial  is  used  to  exactly  constrain  the  transformation.  Given  the  Curtis-Clenshaw  function 
evaluation  points  and  weights  w ,  the  integral  of  a  function  over  an  arbitrarily  dimensioned  hypercube 
is 

Uf{x)dV  =  ^^f(xilj)w^^, 

Vi  7=1 

where  J  is  the  Jacobian  of  the  coordinate  transformation  and  the  weights  can  be  found  from  Burkardt 
(2010).  Because  of  the  non-linear  term  in  the  transformation,  the  Jacobian  is  not  necessarily  constant. 

The  linear  system  for  solution  reconstruction  is  assembled  from  a  stencil  of  Nst  cells.  To  mitigate  the 
effects  of  round-off  error  and  prevent  ill-conditioned  systems,  the  reconstruction  is  computed  on  the 
domain  §-it  =[0,  1]  and  requires  an  additional  linear  transformation  to  go  from  the  fit  domain  to  the 
Curtis-Clenshaw  domain 

^(n)  _  ^cc 

If '  (26) 

where  n  is  the  cell  number  in  the  reconstruction  stencil. 

Equation  (25)  is  written  more  appropriately  using  matrix  notation  as 

(27) 

where  for  the  i-th  cell  in  the  reconstruction  stencil,  w  is  a  column  vector  of  weights  returned  from  the 
Curtis-Clenshaw  quadrature,  P{x-)is  the  coefficient  matrix  for  a  polynomial  fit  of  size  Ncc  x(k+l  f,  and 

C  is  a  column  vector  of  the  unknown  polynomial  coefficients. 

The  complete  linear  system  is  written  as 


det(7(l,J  . 


(25) 


WC  =  V^u,.  (28) 

where  V  is  a  column  vector  of  cell  volumes  and  for  the  ^-th  row  stencil  does 

not  change,  then  this  equation  can  be  solved  prior  to  computing  the  reconstruction  for  computational 
efficiency.  For  all  solution  reconstructions  in  this  work,  polynomial  orders  from  one  to  four  are 
considered. 

k-Exact  Method 

The  k-exact  method  developed  by  Barth  (1990,  1993)  was  designed  to  conserve  the  mean  value  of  the 
cell,  reconstruct  polynomials  of  degree  k  or  less  exactly,  be  compact,  and  be  computationally  efficient. 
For  a  polynomial  of  order  k  there  are  k+1  unknowns  and  Nst  =  k+1.  For  higher  dimensions,  the 
polynomial  is  a  tensor  product  of  one-dimensional  polynomials  with  (k+1)^  unknowns  where  d  is  the 
dimension  and  Nst  =  {k+Vf.  A  centered  stencil  is  used  for  the  interior  of  the  domain  and  a  shifted  stencil 
near  the  computational  boundary.  The  reconstruction  guarantees  that  the  average  value  over  each  cell 
used  in  the  reconstruction  reproduces  the  numerical  solution 
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(29) 


rh  rk 

\=hh^h 

The  k-exact  method  is  solved  in  a  least  squares  sense  by  inereasing  the  size  of  the  steneil  so  that  Equation 
3.25  is  valid  only  for  the  eell  of  interest.  The  method  implemented  here  results  in  the  smallest  possible 
steneil  and  is  equally  sized  in  eaeh  eoordinate  dimension.  The  eonservation  of  the  mean  is  important  in 
the  eontext  of  truneation  error  estimation  so  that 

=  (30) 


ENO  Method 

The  ENO  method  developed  by  Harten  (1987)  was  extended  to  arbitrarily  dimensions  by  Godfrey  et  al. 
(1993).  The  reeonstruetion-via-primitive  ENO  (RP-ENO)  ean  aehieve  arbitrary  high  order  aeeuraey,  but 
is  eomputationally  expensive.  Godfrey  et  al.  (1993)  also  introdueed  a  less  eomputationally  expensive 
dimensionally  split  ENO  (DS-ENO)  method  but  does  not  inelude  eross-derivative  terms.  The  basis  of  the 
ENO  method  is  the  seleetion  of  the  steneil  to  reduee  the  variation  of  the  solution  used  to  eompute  the 
reeonstruetion.  Starting  at  eell  /,  the  differenee  between  the  solutions  to  the  left  (/  -  1)  and  right  (/  +  1)  are 
eompared  and  the  minimum  is  added  to  the  reeonstruetion  steneil.  The  divided  differenee  is  repeated  until 
the  steneil  has  Nst  points. 

RP-ENO  Method 

The  RP-ENO  seheme  eomputes  the  left  and  right  states  of  a  eell  using  progressive  one  dimensional  eurve- 
fits  with  the  adaptive  steneils.  For  a  two-dimensional  reeonstruetion  the  line-averaged  solution  in  the  r|- 
direetion  is  first  eomputed  from  a  one  dimensional  solution  reeonstruetion  where  the  steneil  extends  in  the 
^-direetion 


The  reeonstrueted  solution  is  evaluated  at  the  eell  boundaries  Pi([<^t-i/2^^i+i/2])  in  the  domain. 

Next,  the  solution  is  reeonstrueted  at  the  left  and  right  faees  using  a  one  dimensional  reeonstruetion  of  the 
left  and  right  faee  line  averaged  values  with  a  steneil  that  extends  in  the  rj-direetion.  The  reeonstrueted 
solution  is  then  used  to  eompute  the  flux  at  the  left  and  right  eell  faees.  The  proeess  is  repeated  for  the  top 
and  bottom  faees  of  the  eell  by  eomputing  the  line  averaged  solution  in  the  ^-direetion  first.  A  similar 
proeess  is  followed  for  three  dimensional  reeonstruetions. 

DS-ENO  Method 

The  dimensionally  split  ENO  seheme  follows  the  same  proeedure  exeept  only  one  reeonstruetion  is 
eomputed  for  eaeh  faee  to  eompute  the  line-averaged  solution  in  the  p-direetion  and  the  left  and  right 
faees  and  the  line-averaged  solution  in  the  ^-direetion  at  the  top  and  bottom  faees.  The  average  solution  is 
used  to  eompute  the  flux  instead  of  the  reeonstrueted  loeal  solution  whieh  does  not  eapture  the  eross- 
derivative  terms.  The  eomputational  eost  of  the  DS-ENO  method  is  an  order  of  magnitude  eheaper  than 
the  RP-ENO  and  k-exaet  reeonstruetion  methods;  however,  in  our  ease,  eomputational  expense  is  less  of  a 
eoneern  sinee  the  truneation  error  is  generally  a  one-time  estimate. 


j+l/2 


\\hi^^V)dri. 

i-l/2 


(31) 
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Truncation  Error  Estimation  Results 


Burgers’  Equation 

Burgers’  equation  is  a  useful  sealar  example  problem  to  illustrate  and  initially  test  solution  reeonstruetion 
methods.  Burgers’  equation  is  ehosen  beeause  it  is  analogous  to  the  Navier-Stokes  equations  including  a 
non-linear  convective  term  and  a  linear  diffusive  term.  The  conservative  form  of  Burgers’  equation  is 


(32) 


An  exact  solution  for  a  viscous  shock  is 


u(x)  =  tanh 


(Rq  ] 

[2L  J 


(33) 


where  Re  is  the  Reynolds  number,  L  is  a 
referenee  length  (here  L  =  8),  and  Wre/is  ehosen 
as  two.  The  exaet  solution  with  Re  =  16  is 
shown  in  Figure  1.  This  equation  is  solved 
using  an  explieit,  eell-eentered  finite  volume 
seheme  where 


Figure  1.  Burgers’  Equation  exaet  solution  for  Re  = 
16  and  513  grid  nodes. 


and 


^i+1/2  ^M/2  _  Q 


At 


Ax 


_ul^+uf 
A-+1/2  “  ^ 


(Ar,.,i  +  Ar,.)/2 


(34) 


(35) 


Green’s  theorem  is  used  to  eompute  the  gradient  du/dx  in  Equation  (35)  whieh  is  equivalent  to  a  finite 
differenee  on  a  Cartesian  grid. 
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Figure  2.  Burgers’  Equation  truneation  error 
estimate  for  Re  =  16  and  513  grid  nodes. 


The  full  truneation  error  expression 
eonsists  of  an  infinite  series  and  ean  only  be 
exaetly  ealeulated  using  Equation  (17)  whieh 
requires  the  exaet  solution;  however,  for 
suffieiently  fine  grids,  only  the  leading 
truneation  error  terms  need  to  be  aeeurately 
estimated.  In  general,  the  leading  terms  are  not 
derived  beeause  of  the  eomplexity  of 
diseretization  sehemes,  but  Burgers’  equation 
is  simple  enough  that  this  is  relatively  easy.  On 
a  grid  with  equal  spaeing  the  leading  truneation 
error  terms  are 


t^{u)  =  Ax' 


-u  U  +-UU - vu 

8X  XX  Q  XXX  A  xxxx 

o  Z4 


+  0(Ajc^). 


(36) 


From  Equation  (33)  it  is  elear  that  this  is  a  seeond-order  aeeurate  seheme  and  that  the  next  highest 
truneation  error  terms  are  fourth-order  aeeurate.  For  an  aeeurate  estimate  of  the  leading  terms  the  error 
between  the  estimated  truneation  error  and  the  exaet  truneation  should  deerease  at  at  least  a  fourth-order 
rate. 
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The  discrete  solution  with  Re  =  16 
computed  using  513  grid  nodes  is  reconstructed 
using  the  k-exact  reconstruction  method  with  k 
=  2,  3,  and  4.  The  minimum  value  for  k  is  two 
because  the  highest  derivative  in  the  governing 
equations  is  two  and  any  value  of  k  less  than 
this  would  result  in  a  derivative  term  reducing 
to  zero.  The  truncation  error  is  estimated  from 
the  solution  reconstructions  using  Equation 

(16)  and  the  results  are  shown  in  Figure  2.  The 
estimated  truncation  error  is  compared  to  the 
exact  truncation  error  computed  using  Equation 

(17) .  The  estimated  truncation  error  for  k  =  4  is 
indistinguishable  from  the  exact  truncation 
error  and  the  error  between  the  exact  and 
estimated  truncation  error  decreases  at  a  fourth- 
order  rate.  The  other  two  reconstruction 
methods,  k  =  2  and  3,  do  not  qualitatively 
match  well,  and  the  error  between  the  exact 
and  estimated  truncation  error  decreases  at  only 

a  second-order  rate.  Figure  4  compares  the  maximum  truncation  error  estimate  normalized  by  the 
maximum  exact  truncation  error  for  a  series  of  grids.  The  abscissa  plots  the  grid  spacing  normalized  by 
the  finest  grid  spacing  for  513  nodes.  With  grid  refinement,  it  is  clear  that  the  estimated  truncation  error 
for  k  =  2  and  3  will  never  accurately  represent  the  exact  truncation  error  since  they  asymptote  to  a 
constant  error  with  mesh  refinement.  The  implications  of  this  observation  would  mean  a  minimum 
polynomial  order  k  is  required  to  accurately 
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Figure  4.  Burgers’  equation  normalized  maximum 
truncation  error  estimate  vs.  grid  refinement  for  Re 
=  16. 


estimate  truncation  error  specific  to  each 
discretization  scheme,  and  is  equal  to  the 
highest  derivative  found  in  the  analytic 
truncation  error. 

Euler  and  Navier-Stokes  Equations 

The  same  process  is  applied  for  the  Euler 
and  Navier-Stokes  equations  as  was  done  for 
Burgers’  equation.  The  initial  goal  of  this 
research  is  to  determine  the  minimum  order 
polynomial  reconstruction  required  to 
accurately  estimate  the  truncation  error  for  the 
Euler  and  Navier-Stokes  equations.  Solutions 
and  truncation  error  estimates  are  computed 
using  three  different  grid  families  with  the 
coarsest  grid  being  17x17  nodes  and  the  finest 
grid  being  257x257  nodes.  An  example  of  each 
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Figure  3.  Curvilinear  grid,  33x33  grid  nodes, 
is  shown  in  Figure  3.  The  Cartesian  grid  is  the  simplest  and  is  expected  to  be  the  easiest  to  estimate 
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truncation  error.  The  complexity  is  increased  on  the  skewed  grid;  however,  the  linear  grid  transformation 
can  exactly  represent  the  grid  distribution.  The  final  grid  is  a  curvilinear  grid  with  a  sinusoidal 


K=2 


Figure  5.  Euler  energy  truncation  error  estimation  comparisons  on  the  curvilinear  mesh  for  the 
supersonic  solution  with  a  viscosity  of  p  =  IPa-s. 


distribution  which  cannot  be  exactly  represented  by  the  linear  grid  transformation. 

Truncation  error  estimates  for  Euler  and  Navier-Stokes  solutions  are  computed  using  the  curvilinear 
grid.  The  different  methods  include  k-exact,  RP-ENO,  DS-ENO,  the  coarse  grid  method  used  by  Fulton 
(2003),  with  and  without  the  correction  term,  and  the  fine  grid  method  used  by  Venditti  and  Darmofal 
(2003),  with  and  without  the  correction  term,  for  reconstructions  with  k  =  [1,2,3, 4].  The  effectivity  index 
(Ainsworth  and  Oden,  2000)  is  used  to  evaluate  the  truncation  error  estimates  which  is  the  truncation 

error  estimate  normalized  by  the  exact  truncation  error,  9^2  ~  lk/z(^/z^/z)|L2 ‘  ^  shows  the 

effectivity  index  for  the  energy  truncation  error  for  the  supersonic  Euler  solution,  and  Figure  6  shows  the 
effectivity  index  for  the  energy  truncation  error  for  the  subsonic  Navier-Stokes  solution.  The  subsonic 
Euler  and  Supersonic  Navier-Stokes  results  are  very  similar  to  the  subsonic  Navier-Stokes  and  supersonic 
Euler  results,  respectively.  The  energy  equation  was  shown  because  it  is  the  most  complex  truncation 
error  term;  however,  the  truncation  error  for  the  other  equations  behave  in  a  similar  manner. 
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Figure  6.  Navier-Stokes  energy  truneation  error  estimation  eomparisons  on  the  eurvilinear  mesh  for 

the  subsonie  solution  with  a  viseosity  of  |i  =  IPa-s. 


For  both  the  Euler  and  Navier-Stokes  equations,  the  k-exaet  reeonstruetion  method  is  generally  the 
best  performing  method.  For  k=  1,  the  truneation  error  estimates  deerease  at  a  first-order  rate  where  the 
exaet  truneation  error  deereases  at  a  seeond-order  rate.  The  lower  rate  results  in  the  estimated  truneation 
error  diverging  from  the  exaet  truneation  error  with  grid  refinement.  For  k  =  2  and  greater,  the  k-exaet 
truneation  error  estimates  are  very  aeeurate  and  outperform  the  other  methods  with  only  a  few  exeeptions. 
Based  on  our  previous  results  with  Burgers’  equation,  it  was  expeeted  that  a  higher-order  reeonstruetion 
would  be  required  for  the  Navier-Stokes  equations;  however,  all  results  show  that  k  =  2  is  suffieient  for 
both  the  Euler  and  Navier-Stokes  equations.  A  possible  explanation  for  this  result  is  due  to  the 
diseretization  seheme  for  the  diffusion  terms.  The  derivatives  are  eomputed  using  Green’s  theorem  whieh 
reduees  to  a  seeond-order  eentral  differenee  for  one  dimensional  solutions  on  an  evenly  spaeed  grid.  The 
truneated  terms  for  the  seeond-order  eentral  differenee  is  0(Ax^)  and  (9(Av'^).  The  eonveetive  terms  have 
higher  order  terms  whieh  are  0(Ax^)  and  0(Ax^).  This  would  mean  that  the  eonveetive  truneation  error 
terms  will  dominate  the  error  in  the  truneation  error  estimate.  A  higher  viseosity  Navier-Stokes  solution 
was  eomputed  to  try  to  inerease  the  dominanee  of  the  diffusive  terms.  The  manufaetured  solution  souree 
terms  were  eompared,  and  the  diffusive  terms  were  on  the  same  order  of  magnitude  as  the  eonveetive 
terms.  Further  evidenee  that  k  =  2  is  suffieient  is  shown  by  the  divergenee  of  the  ENO  methods  for  the 
supersonie  solutions  shown  in  Figure  7  and  Figure  8  whieh  are  missing  derivatives  due  to  the  lower 
dimension  reeonstruetion.  If  the  diffusive  terms  were  insignifieant,  the  ENO  methods  would  be  aeeurate 
as  shown  in  Figure  7.  It  might  be  possible  that  these  terms  eould  dominate  for  highly  diffusive  areas  of  a 
flow  requiring  a  minimum  of  k  =  3;  however,  there  is  no  evidenee  to  suggest  that  k  =  2  is  not  suffieient. 

For  the  Euler  equations  and  k=  1,  the  ENO  truneation  error  estimation  methods  deerease  at  a  first- 
order  rate  similar  to  the  k-exaet  method  and  results  in  an  effeetivity  index  that  diverges  with  grid 
refinement  (not  shown  in  the  figure).  For  k  =  2,  the  ENO  truneation  error  estimation  methods  deerease  at 
a  third-order  rate  resulting  in  a  substantial  under  predietion  of  the  truneation  error.  The  exaet  reason  for 
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Figure  7.  Navier-Stokes  energy  truneation  error  estimation  eomparisons  on  the  eurvilinear  mesh  for 
the  subsonie  solution  with  a  viseosity  of  |i  =  IPa-s. 

this  behavior  is  not  yet  understood;  however,  the  ENO  sehemes  are  used  for  solution  reeonstruetion  where 
k  =  2  is  third  order  aeeurate.  The  truneation  error  estimate  for  k  =  2  eould  be  an  estimate  of  the  truneation 
error  not  for  the  eurrent  MUSCL  extrapolation  upwind  method  but  an  approximation  of  the  third-order 
aeeurate  ENO  seheme;  however,  the  truneation  error  estimates  for  k  =  3  and  k  =  4  results  in  aeeurate 
truneation  error  estimates  with  only  slightly  higher  error  than  the  k-exaet  method.  For  the  Navier-Stokes 
equations,  the  subsonie  solution  follows  a  similar  trend  as  the  Euler  equations;  however,  for  the 
supersonie  solution  the  ENO  methods  do  not  aeeurately  estimate  the  truneation  error.  The  estimate  is 
zeroth  order  aeeurate  beeause  the  ENO  methods  eannot  represent  the  eross-derivative  terms.  The 
differenee  between  the  subsonie  and  supersonie  Navier-Stokes  solution  (shown  in  Figure  7)  is  thought  to 
be  due  to  the  relative  magnitude  of  the  eonveetive  and  diffusive  terms.  The  supersonie  solution  has  mueh 
larger  eonveetive  terms  due  to  the  larger  magnitude  veloeity  and  veloeity  gradient  and  the  zeroth  order 
terms  are  more  apparent.  The  ENO  sehemes  ean  be  used  for  truneation  error  estimation;  however, 
aeeuraey  would  suffer  in  diffusive  dominated  flow  regimes  sueh  as  boundary  layers  or  shear  flows.  The 
effeets  of  a  more  dominant  diffusive  term  is  shown  in  Figure  8  eompared  to  Figure  7.  (Due  to  stability 
issues  with  the  higher  viseosity,  only  three  grid  levels  eonverged  iteratively.)  The  results  are  similar, 
exeept  the  truneation  error  estimated  from  the  solution  eomputed  with  a  viseosity  of  50Pa-s  diverges  more 
quiekly  than  the  truneation  error  estimated  from  the  solution  eomputed  with  a  viseosity  of  IPa-s. 

For  the  eoarse  grid  truneation  error  estimation  methods,  solution  reeonstruetion  is  used  to  interpolate 
the  eoarse  grid  truneation  error  estimate  baek  to  the  eomputational  mesh.  The  uneorreeted  eoarse  grid 

truneation  error  estimation  method  is  off  by  a  faetor  of  l/(r^^  -1) ;  however,  the  eorreeted  method 
aeeurately  estimates  the  truneation  error.  The  most  aeeurate  reeonstruetion  seheme  for  the  eorreeted 
method  uses  a  k-exaet  reeonstruetion  method  with  k  =  1  for  the  Euler  equations.  For  the  Navier-Stokes 
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equations,  k=\  underestimates  the  error  for  the  subsonie  solution  and  the  more  diffusive  supersonie 
solution;  therefore,  for  the  Navier-Stokes  equations  k  =  2  is  reeommended. 


Figure  8.  Navier-Stokes  energy  truneation  error  estimation  eomparisons  on  the  eurvilinear  mesh  for 
the  supersonie  solution  with  a  viseosity  of  p  =  50Pa-s. 


For  the  fine  grid  truneation  error  estimation  methods,  solution  reeonstruetion  is  used  to  interpolate  the 
eomputational  solution  onto  a  finer  mesh  that  is  refined  by  a  faetor  of  two  in  all  eoordinate  direetions.  The 
uneorreeted  fine  grid  truneation  error  estimation  method  underestimates  the  truneation  error.  Correeted  by 

a  faetor  of  l(r^^  -1)  results  in  an  aeeurate  truneation  error  for  k  =  2  or  higher  for  both  the  Euler 
equations  and  the  Navier-Stokes  equations.  The  fine  grid  method  with  the  eorreetion  is  one  of  the  most 
aeeurate  truneation  error  estimators  evaluated.  The  method  is  more  aeeurate  for  k  =  2  than  the  ENO  and 
eoarse  grid  methods  and  indistinguishable  from  the  k-exaet  method  for  all  solutions  exeept  the  highly 
diffusive  Navier-Stokes  solutions  in  whieh  the  k-exaet  truneation  error  estimate  is  more  aeeurate.  It  is 
important  to  note  that  the  test  problems  are  very  smooth.  The  eorreetion  term  assumes  that  the  truneation 
error  deereases  at  the  formal  order  of  aeeuraey  whieh  may  not  be  the  ease  for  highly  non-asymptotie 
solutions,  solutions  with  poor  grid  quality,  or  solutions  with  singularities/diseontinuities. 

Contour  plots  of  the  Navier-Stokes  energy  truneation  error  estimates  for  the  supersonie  solution  on 
the  33x33  grid  are  shown  in  Figure  9.  The  reeonstruetion  order  for  eaeh  method  was  ehosen  from  the  best 
results  from  the  previous  results  shown.  All  eontour  plots  have  the  same  eontour  levels.  All  eontour  plots 
eompare  qualitatively  well  to  the  exaet  truneation  error.  The  k-exaet  and  fine  grid  method  with  the 
eorreetion  term  eompare  very  well  to  the  exaet  truneation  error.  The  eoarse  grid  method  with  the 
eorreetion  faetor  also  eompares  well  but  has  a  few  peaks  that  are  not  present  in  the  exaet  truneation  error. 
The  ENO  methods  qualitatively  mateh  well,  but  the  truneation  error  estimates  are  not  as  smooth.  The 
truneation  error  for  the  ENO  methods  on  the  33x33  grid  is  the  most  aeeurate  estimate  for  the  series  of 
grid  levels  (see  Figure  l,k=  3).  On  the  65x65  grid  the  truneation  error  estimate  begins  to  diverge. 
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Figure  9.  Truncation  error  estimate  contours  for  Navier-Stokes  energy  equation  on  a  33x33 
curvilinear  grid. 


In  the  presence  of  shocks,  the  order  of  accuracy  of  the  numerical  scheme  reduces  to  first-order 
accurate.  The  same  results  were  computed  for  first-order  Euler  and  Navier-Stokes  solutions  to  determine 
the  minimum  required  reconstruction  order  for  the  truncation  error  estimation  methods.  The  results  are 
shown  in  Figure  10.  The  results  show  that  the  reconstruction  order  is  one  less  than  what  is  required  for 
second-order  accurate  solutions  for  both  the  Euler  and  Navier-Stokes  equations.  All  methods  are  accurate 
for  k  =  1  except  the  ENO  schemes  which  require  k  =  2.  The  coarse  grid  correction  term  reduces  to  one 
(i.e.,  it  is  not  needed)  for  first-order  solutions  and  the  fine  grid  correction  term  is  two  as  expected. 
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Figure  10.  Navier-Stokes  energy  truneation  error  estimation  eomparisons  on  the  eurvilinear  mesh  for 
the  supersonie  first-order  solution. 


The  results  for  all  the  truneation  error  estimation  methods  are  summarized  as  follows.  The  two  best 
performing  methods  are  the  k-exaet  and  the  fine  grid  method  with  the  eorreetion  term.  The  k-exaet 
method  does  not  require  the  assumption  of  an  order  of  aeeuraey  whieh  is  a  signifieant  advantage, 
espeeially  for  more  praetieal  applieations  where  the  formal  order  of  aeeuraey  may  not  be  aehieved.  For 
the  first-order  aeeurate  methods,  the  ENO  methods  minimum  required  reeonstruetion  order  for  the 
Navier-Stokes  equations  showed  k  =  2  is  suffieient.  This  is  beeause  the  diffusion  terms  (evaluated  using 
Green’s  theorem)  are  still  seeond-order  aeeurate  while  the  eonveetive  terms  are  first-order  aeeurate.  The 
eonveetive  error  terms  dominate  the  truneation  error  and  results  in  truneation  error  estimates  that  are 
nearly  identieal  to  the  Euler  equations.  This  result  suggest  that  the  ENO  methods  eould  be  used  in  the 
vieinity  of  flow  singularities  sueh  as  shoek  or  eontaet  diseontinuities,  where  the  order  of  aeeuraey  reduees 
to  first  order  or  below  (Banks  et  al.,  2008). 

There  are  several  uses  for  truneation  error  whieh  inelude  diseretization  error  estimation  and  higher 
order  solution  eorreetion.  The  possible  diseretization  error  estimation  methods  inelude  the  error  transport 
equations  (Zhang  et  al.,  2002;  Qin  and  Shih,  2002),  defeet  eorreetion  (Pereyra,  1965;  Stetter,  1978;  Skeel, 
1986),  and  adjoint  methods  (Giles  and  Pieree,  2000).  Defeet  eorreetion  is  the  least  eode  intrusive  to 

implement  as  the  estimated  truneation  error  is  added  as  a  souree  term,  L^(u)  =  .  The  solution  u 

is  an  estimate  of  the  exaet  solution,  and  the  diseretization  error  is  estimated  by  -u^—u  .  An  example 

of  defeet  eorreetion  is  shown  in  Figure  1 1  eomputed  using  all  truneation  error  methods  for  the  supersonie 
Euler  solution  on  the  eurvilinear  grid.  A  sliee  is  taken  through  the  eenter  of  the  domain  at  the  j  =  16  eell 
index.  The  k-exaet  and  ENO  truneation  error  estimation  methods  are  the  most  aeeurate  with  the  eorreeted 
eoarse  grid  method  performing  well.  The  uneorreeted  eoarse  grid  method  is  not  aeeurate  and 
overestimates  the  error  signifieantly.  The  fine  grid  methods  do  not  eapture  the  shape  of  the  diseretization 
error  very  well. 
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Figure  11.  Defect  correction  example  using  the  supersonic  Euler  solution  on  the  33x33  grid. 


The  second  solution  w  is  a  higher  order  accurate  solution  and  can  be  used  as  a  solution  instead  of 
estimating  the  discretization  error.  Again,  for  the  supersonic  Euler  equations,  the  k-exact  method  is  used 
to  compute  higher-order  solutions  shown  in  Figure  12.  The  expected  order  of  accuracy  is  third-order 
because  the  k-exact  method  is  estimating  the  leading  truncation  error  terms  which  are  O(Ax^).  The  order 
of  accuracy  begins  very  high  and  appears  to  be  approaching  third-order  for  the  finest  grid.  The  higher- 
order  discretization  error  is  not  as  asymptotic  as  the  second-order  solution;  however,  the  discretization 
error  on  the  finest  grid  is  on  the  order  of  1x10'^  %  error.  The  discretization  error  in  the  corrected  solution 
is  lower  for  all  solution  variables  except  for  the  pressure  which  is  higher  than  the  second-order  solution, 
but  because  of  the  higher  order  of  accuracy  quickly  decreases.  The  solver  does  not  have  any  higher-order 
capability,  the  higher-order  results  come  from  the  truncation  error  estimate  only. 
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Figure  12.  Discretization  error  in  the  defect  correction  solution  discretization  error  (left)  and  the  order 
of  accuracy  (right)  computed  using  the  k-exact  truncation  error  estimation  method  for  the  supersonic 
Euler  solution. 


VII.  Residual-Based  Discretization  Error  Estimation 

In  general,  the  focus  of  the  discretization  error  estimation  is  local  solution  estimation.  The  local 
solution  variables  for  the  compressible  Euler  and  Navier-Stokes  equations  are  density,  pressure,  and 
velocity  components.  To  evaluate  the  accuracy,  contour  plots  are  used;  however,  for  a  more  quantitative 
analysis  discrete  Li-  and  L2-norms  are  used.  The  Li-norm  is  used  when  shocks  are  present  in  the  solution 


and  the  L2-norm  is  used  otherwise 


Z 


'Afti 


(37) 


(38) 


To  compare  the  estimated  errors  to  the  exact  errors,  the  effectivity  index  is  used  (Ainsworth  and  Oden, 

2000) 


^L2  ~  ■ 


\\L2 


\\L2 


(39) 


where  6*^  is  an  estimated  error  and  6*^  is  the  exact  error.  If  the  estimated  error  is  accurate  0l2  «  1.  For 

asymptotically  accurate  error  estimates  the  effectivity  index  should  approach  one  as  the  computational 
grid  is  refined.  The  effectivity  index  is  used  to  evaluate  local  truncation  error  and  local  discretization  error 
estimates  for  Li-  and  L2-norms  where  appropriate. 

The  residual -based  discretization  error  estimation  methods  that  we  will  investigate  include  error 
transport  equations,  defect  correction,  and  adjoint  methods.  All  three  of  these  methods  have  both  a 
continuous  and  a  discrete  implementation  and  use  discrete  or  continuous  residuals  that  can  be  related  back 
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to  the  truncation  error.  The  error  transport  equations  and  defect  correction  both  involve  solutions  of 
additional  governing  equations  on  the  same  mesh  with  residuals/truncation  errors  acting  as  local  source 
terms;  however,  the  error  transport  equations  involve  an  additional  linearization  step  while  defect 
correction  relies  on  the  “nearness”  of  a  nearby  problem  to  the  original  one.  Here  we  reformulate  all  three 
residual-based  error  estimations  methods  within  the  new  GTEE  framework  of  Equation  (10). 

Error  Transport  Equations 

Using  the  new  GTEE  framework  given  by  Equation  (10),  the  continuous  discretization  error  transport 
equation  can  be  found  as  follows.  First,  insert  into  the  GTEE  and  subtract  the  restriction  of  the  PDE 

=  0)  to  give 


I%I,u,)-I>'L{u)  =  -TE,{I,u,)  (40) 

where  we  have  assumed  that  =  L{u^)  =  0  (i.e.,  that  the  restriction  of  the  prolongation 

returns  the  discrete  value).  Prolonging  Equation  (40)  to  a  continuous  space  and  making  the  approximation 
»  L(w)  (which  is  not  generally  true,  but  may  be  a  reasonable  approximation)  results  in: 


(41) 


If  the  equations  are  linear,  or  if  we  employ  simple  linearization,  then  L{I^Uj^)-  L{u)  =  L{I^u^  -u)  =  L{s) 
thus  resulting  in  the  continuous  discretization  error  transport  equation 


L(e)  =  -/,7E, (/,»,) 


(42) 


which  for  our  ID  Burgers’  equation  example  becomes: 

^ds  d^£  j  rrj,  X 

^^-^^  =  -hTEh{h^h)- 
ax  ax 

Here  we  propose  a  more  advanced  linearization  by  noting  that  the  operators  for  Burgers’  equation  in 
Equation  (41)  can  be  written  out  as 


E{hH)-E{u)  =  IhUh 


dih^h)  ^du  ^d^u 

V  2  ^  ^  2  ■ 

dx  dx  dx  dx 


(43) 


Using  the  definition  of  the  discretization  error  from  Equation  (1 1)  to  give  =  u+£  and  substituting 
into  Equation  (43)  gives 


L{I -  L(u)  =  (u  +  s) 


d{Ii^Ui^)  ^du  d^u 


dx 


dx^ 


dx  dx^ 


or  simply 
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(44) 


ax  ax 


ds  dh  d(I,u,)_-  dil,u,) 


=  L{e)  +  s- 


dx  dx 

where  the  overbar  on  the  PDE  operator  indieates  a  linearization.  With  this  advaneed  linearization,  the 
eontinuous  diseretization  error  transport  equation  for  Burgers’  equation  beeomes 

^  —--1  TF  (1  n  \ 


^ds 

u - V- 

dx  dx 


dx 


It  is  important  to  inelude  the  additional  souree  term  sinee  both  the  truneation  error  and  the  diseretization 
error  s  will  have  the  same  formal  order  of  aeeuraey.  As  will  be  shown  later,  omitting  this  term  results  in 
error  estimates  whieh  are  ineorreet  by  a  fixed  faetor  (e.g.,  they  are  always  30%  low  regardless  of  the  level 
of  mesh  refinement).  Shih  and  Williams  (2009)  address  this  linearization  issue  indireetly  by  formulating 
the  residual  on  a  finer  mesh  rather  than  employing  the  true  truneation  error. 

A  similar  diserete  diseretization  error  transport  equation  ean  be  found  by  inserting  u  into  the  GTEE 
and  subtraeting  the  diserete  equations  =  0)  to  give 


(45) 


If  the  equations  are  linear,  or  if  we  employ  a  simple  linearization,  then 

Lh {u^ )-Lj^ (I^u)  =  -  I^u)  =  )  thus  resulting  in  the  diserete  diseretization  error  transport 

equation 


(46) 


whieh  for  the  simple  2^^  order  aeeurate  finite  differenee  diseretization  of  ID  Burgers’  equation  beeomes: 


_(+! - L_L_y_^±i - L - L±  -  -TEAu)  . 

'  2Ax  Ax^  ' 


(47) 


The  more  advaneed  linearization  has  (I  u)  =  )  +  S- 


2Ax 


resulting  in 


(48) 


2 Ax  Ajc  ""  "  '  2Ajc 

Note  that  the  presenee  of  the  exaet  solution  to  the  PDEs  in  the  above  error  transport  equations  ean  always 
be  approximated  by  the  numerieal  solution  and  the  error  estimate:  u.-I^u-u.-  6*. . 

Sinee  we  already  have  an  implieit  solution  eapability  in  our  eode,  we  have  found  that  the  following 
alternative  linearization  provides  an  effieient  form  of  the  error  transport  equations.  The  diserete  error 
transport  equations  are  derived  by  substituting  the  exaet  solution  to  the  eontinuous  governing  equation 
into  the  GTEE  given  in  (10),  subtraeting,  L^(w^)  =  0  ,  and  noting  that  L{u)  =  0 : 


LM-L,{l‘'d)  =  -T,{u) 


(49) 
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The  second  term  on  the  left-hand  side  can  be  linearized  in  a  Taylor  series  expanded  about  the  numerical 
solution  u^\ 


(50) 


Inserting  this  equation  into  Equation  (49),  we  obtain  the  final  error  transport  equation  result: 


(51) 


du 


The  derivative  is  the  global  Jacobian  of  the  linearized  discrete  operator  which  is  required  to  implement  an 


implicit  scheme.  The  Jacobian  is  commonly  available  in  flow  solvers;  however,  the  Jacobian  required  to 
solve  the  error  transport  equations  should  be  the  full  higher  order  (e.g.,  second-order)  left-hand  side  where 
only  the  first-order  left-hand  side  may  be  required  for  an  implicit  solve.  The  above  form  of  the  error 


transport  equation  is  approximate  so  the  exact  discretization  error  may  not  result  if  the  exact  truncation 
error  is  used;  however,  the  error  in  Equation  (51)  decreases  at  a  much  higher  rate  and  would  introduce 
non-negligible  error  only  in  the  case  of  very  large  discretization  errors. 

Defect  Correction 

Defect  correction  methods  were  originally  developed  nearly  50  years  ago  in  order  to  improve  the 
accuracy  of  (or  estimate  the  numerical  error  in)  numerical  solutions  to  ordinary  differential  equations 
(e.g.,  Zadunaisky,  1966;  Stetter,  1978).  More  recently,  defect  correction  has  been  used  in  CFD  to  increase 
the  accuracy  of  order  spatial  discretization  schemes  (Layton  et  ah,  2002;  Naumovich  et  ah,  2010).  We 
are  interested  in  using  defect  correction  methods  to  estimate  the  discretization  error  in  CFD  solutions. 
There  are  two  main  types  of  defect  correction  (Skeel,  1986):  continuous  defect  correction  (also  called 
differential  correction)  and  discrete  defect  correction  (also  called  difference  correction).  We  will  present 
the  main  ideas  behind  our  approach  to  continuous  defect  correction,  then  simply  explain  how  discrete 
defect  correction  would  differ. 

Defect  correction  corrects  the  original  numerical  solution  by  re -solving  the  original  problem  with  the 
truncation  error  estimate  added  to  the  right-hand  side.  Defect  correction  requires  minimal  code  intrusion. 
The  exact  defect  correction  problem  is  defined  as 


L^{l‘'u)  =  Th{u) 


(52) 


If  the  exact  truncation  error  T^{u)  is  added  as  a  source  term,  the  resulting  numerical  solution  is  the  exact 

solution  to  the  PDEs  or  integral  equations  restricted  to  the  computational  grid.  Our  defect  correction 
implementation  takes  the  form 


(53) 


where  u  and  the  discretization  error  estimate  is 


(54) 
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Defect  correction  requires  a  second  numerical  solution  on  the  original  computational  mesh;  however, 
computational  expense  is  reduced  because  the  original  numerical  solution  to  the  primal  problem  is  used  to 
initialize  the  second  simulation,  thus  providing  a  very  good  starting  solution. 

Adjoint  Methods 

Adjoint  methods  provide  a  means  by  which  discretization  error  estimation  and  mesh  adaptation  can 
be  targeted  to  a  solution  functional  (e.g.,  lift  or  drag).  We  begin  our  discussion  of  adjoint  methods  by 
reinterpreting  adjoint  methods  in  the  GTEE  framework.  Consider  a  scalar  solution  functional  J{u)  and 

its  discrete  counterpart  .  We  may  be  interested  in  estimating  the  discretization  error  in  this 

functional 


(55) 


assessing  the  sensitivity  of  this  functional  to  local  values  of  the  truncation  error,  or  both.  For  the 
continuous  adjoint,  we  first  form  the  Lagrangian  using  the  inner  product  (e.g.,  (/,  ^)  =  j  / {x)g{x)dx ) 


f(w,^)  =  /(w)-{%  L{u)) 

Next,  linearize  both  J  and  L  about  the  general  function  u: 


(56) 


ou 


_  ^  dL 

(u  -u)  +  HOT  and  L(m)  =  L(m)  h - 

du 


{u-u)  +  HOI . 


Inserting  these  linearizations  into  the  right-hand  side  of  Equation  (56)  and  neglecting  the  higher-order 
terms  gives 


f(M,^)  =  y(M)-{^,L(M))  + 


dJ 

du 

u  \ 

./J 

(u-u) . 


(57) 


Replacing  the  general  function  u  in  Equation  (57)  with  the  prolongation  of  a  numerical  solution,  i.e., 
u  =  have 


where  the  adjoint  (dual)  problem  is  defined  as 


dJ 


du 


-  'F. 


du 


4% , 


(58) 


dJ 

=  { 

) 

du 

h^h 

1  du 

h^h  1 

(59) 


Once  the  adjoint  problem  is  solved  for  the  term  in  brackets  in  Equation  (58)  will  be  equal  to  zero.  By 
noting  that  the  PDE  is  equal  to  zero  over  the  entire  domain  (i.e.,  L{u)  =  0 ),  we  can  use  Equation  (56)  to 
replace  the  Lagrangian  with  =  J{u)  in  Equation  (58)  to  give 
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(60) 


where  L{I^u^)  is  simply  the  eontinuous  residual  (an  estimate  of  the  truneation  error).  If  the  funetional  J 
involves  an  integral,  then  the  numerieal  integration  errors  ean  be  estimated  as 

(61) 

Combining  Equations  (60)  and  (61)  with  the  definition  of  the  diseretization  error  in  the  funetional  given 
by  Equation  (55),  we  have 

=  +  =  ^integ  +  Lilht^h  )) 


or  simply 


~  J{u)  -  •  (62) 

The  diseretization  error  in  the  solution  funetional  is  thus  due  to  1)  the  integration  error  (whieh  ean  be 
redueed  simply  by  improving  the  numerieal  quadrature  used  in  evaluating  the  integral)  and  2)  the  loeal 
truneation  error  (or  residual)  weighted  by  the  adjoint  sensitivities  W.  The  solution  to  the  adjoint  problem 
thus  provides  the  sensitivity  of  the  solution  funetional  to  loeal  truneation  error  sourees.  A  diserete  adjoint 
eould  be  derived  similarly  whieh  uses  the  diserete  residual: 

e,=JM-J(u)  =  -  'P  "i.  (/‘ff)  (63) 

where  is  a  row  veetor  eontaining  the  adjoint  sensitivities  and  L^{I^u)  is  a  eolumn  veetor  eontaining 

the  loeal  diserete  residual  (i.e.,  the  truneation  error)  for  eaeh  eell  in  the  mesh.  Both  the  eontinuous  and 
diserete  forms  of  the  adjoint  method  will  be  investigated  within  the  GTEE  framework.  Methods  for 
formulating  the  sensitivities  needed  for  the  adjoint  methods  will  be  examined  ineluding  analytieal 
derivative  evaluation,  numerieal  derivative  evaluation,  and  eomplex  numbers  (Squire  and  Trapp,  1998). 

VIIL  Local  Error  Estimation  Results 

Local  Error  Estimates:  NACA  0012  Airfoil 

For  the  loeal  error  estimates,  the  truneation  error  is  estimated  using  the  k-exaet,  least  squares  (ESQ), 
DSENO  single  grid  methods  as  well  as  the  eoarse  grid  method  (with  and  without  the  eorreetion  term)  and 
the  fine  grid  method  (with  and  without  the  eorreetion  term).  Defeet  eorreetion  and  ETEs  are  solved  with 
eaeh  of  the  truneation  error  estimates. 
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The  truncation  error  estimates  for  the  subsonic  airfoil  at  zero  degrees  angle  of  attack  are  shown  in 
Figure  14.  The  k-exact  and  fine  grid  method  with  the  correction  term  capture  the  truncation  error  around 
the  leading  edge  of  the  airfoil.  The  least  squares  method  compares  well  also  with  the  exception  of  a  large 
truncation  error  estimate  at  the  leading  edge.  The  DS-ENO  method  results  in  a  very  noisy  truncation  error 
estimate.  The  discretization  error  in  pressure  is  shown  in  Figure  13.  The  accuracy  of  the  truncation  error 
estimates  directly  correspond  to  the  accuracy  of  the  discretization  error  estimate.  The  k-exact  estimate 
compares  well  to  the  exact  discretization  error  and  the  DS-ENO  method  is  relatively  noisy. 
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Figure  14.  Estimated  truncation  error  for  the  mass  eqn.  for  the  NACA  0012  airfoil  at  M=0.5  and  a=0. 


Exact 


K-exact 


LSQ 


DS-ENO 


Richardson  Extrapolation 


Coarse  Corrected 


Fine  Corrected 


0.5 

0^ 


I  452.554 
;  269.103 
^  85.6526 
.  -97.798  0.5 

‘-281.249  ^ 


Ilk.. 


05 

X 


05 

X 


_  1 

■  452.554 

452.554 

269.103 

269.103 

'  85.6526 

L.  , 

85.6526 

t  -97.798  0.5 

-97.798 

*“-281.249 

-281  249 

0 

V* 

p. 

i 

6  0.5 

1 

X 

Figure  13.  Estimated  defect  correction  discretization  error  estimates  for  pressure  for  the  NACA  0012 
airfoil  at  M=0.5  and  a=0. 
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To  compare  the  difference  between  defect  correction  and  the  ETEs,  the  exact,  defect  correction  estimate, 
and  two  EXE  estimates  are  shown  in  Figure  15.  One  EXE  estimate  is  computed  using  a  Jacobian 
computed  for  first-order  discrete  equations  and  the  other  is  computed  using  a  Jacobian  computed  for  the 
second  order  discrete  equations.  The  results  are  very  nearly  identical  and  only  differ  by  fractions  of  a 
percent.  The  pressure  on  the  upper  wall  is  also  compared  in  Figure  16.  The  results  are  nearly  identical  for 
the  EXE  estimates.  The  result  is  fortuitous  because  a  first  order  Jacobian  is  often  readily  available  when 
using  an  implicit  solver.  The  solution  time  for  the  ETEs  using  the  first  order  Jacobian  is  negligible 
compared  to  the  expense  of  the  defect  correction  estimate.  For  comparison,  the  primal  solution  took  24 
minutes  while  using  the  k-exact  truncation  error  estimation  method,  the  defect  correction  error  estimate 
requires  about  nine  minutes  to  solve  while  the  first-order  EXE  took  less  than  a  second  to  compute  with 
three  iterations.  The  second-order  EXE  cost  about  five  seconds  to  compute  with  seven  iterations.  The 
first-order  EXE  is  so  inexpensive  because  the  equations  are  linear.  The  difference  between  the  first-  and 
second-order  EXE  is  negligible  for  all  applications  so  only  the  first-order  EXE  results  are  shown.  The 
EXE  estimates  slightly  underestimate  the  discretization  error  at  the  leading  edge  compared  to  defect 
correction.  Both  residual-based  methods  accurately  estimate  the  discretization  error  while  Richardson 
extrapolation  misses  the  sign  and  magnitude  near  the  leading  edge  significantly. 


Exact  Defect  Correction 


ETE  (1  St  order  Jacobian)  ETE  (2nd  Order  Jacobian) 


Figure  15.  Estimated  defect  correction  and  ETE  discretization  error  estimates  for  pressure  for  the 
NACA  0012  airfoil  at  M=0.5  and  a=0  using  the  k-exact  reconstruction  method  (solid  is  the 
benchmark  solution). 
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Discretization  Error  Comparison 


Figure  16.  Estimated  defeet  eorreetion  and  EXE  diseretization  error  estimates  for  pressure  along  the 
upper  side  of  the  the  NACA  0012  airfoil  at  M=0.5  and  a=0  using  the  k-exaet  reeonstruetion  method 
(solid  is  the  benehmark  solution). 


Figure  17  and  Figure  18  eompare  the  pressure  on  the  wall  for  defeet  eorreetion  and  EXE 
diseretization  error  estimates  eompared  to  the  exaet  diseretization  error.  For  all  eases,  the  EXE  estimate 
eompares  well  to  the  defeet  eorreetion  estimate.  Xhere  are  slight  differenees  near  diseretization  error 
peaks,  and  overall,  the  defeet  eorreetion  estimate  is  more  aeeurate  beeause  the  estimate  results  from 
nonlinear  equations  eompared  to  the  linearized  EXEs.  In  terms  of  truneation  error  estimation,  the  k-exaet 
method  is  the  most  aeeurate  with  only  a  minor  differenee  in  peak  error  near  the  leading  edge.  Xhe  ESQ, 
DS-ENO,  and  fine  grid  method  with  the  eorreetion  faetor  eompare  well  for  the  entire  airfoil  with  larger 
differenees  near  the  leading  edge  of  the  airfoil.  All  methods  eompare  better  than  Riehardson  extrapolation 
exeept  for  the  eoarse  grid  methods  whieh  are  inaeeurate  over  the  entire  airfoil. 
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Figure  17.  Estimated  defect  correction  discretization  error  estimates  for  pressure  along  the  upper  side 
of  the  NACA  0012  airfoil  at  M=0.5  and  a=0  (solid  is  the  benchmark  solution). 
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Figure  18.  Estimated  ETE  discretization  error  estimates  for  pressure  along  the  upper  side  of  the 
NACA  0012  airfoil  at  M=0.5  and  a=0  (solid  is  the  benchmark  solution). 


IX.  Functional  Error  Estimation  Results 

As  alternatives  to  the  adjoint  method,  whieh  ean  only  provide  a  DE  for  a  single  fiinetional  output  per 
dual  solve,  both  defeet  eorreetion  methods  (DC)  (Pereyra,  1965,  1967,  1968)  and  the  error  transport 
equations  (ETEs)  (Phillips  and  Roy,  2011)  provide  loeal  DE  estimates  driven  by  TE  estimates  and  ean  be 
used  to  estimate  error  in  any  funetional  of  interest  for  a  single  solve.  DC  methods  treat  the  TE  estimate  as 
a  souree  term  to  drive  the  primal  solver  towards  a  higher  order  solution.  This  TE  souree  term  ean  be 
thought  of  in  a  similar  manner  to  the  souree  term  utilized  by  the  method  of  manufaetured  solutions 
(MMS)  (Roaehe,  2009;  Oberkampf  and  Roy,  2010)  to  drive  the  diserete  governing  equations  towards  a 
preseleeted  solution.  If  the  exaet  TE  is  known,  then  the  diserete  governing  equations  ean  be  driven 
towards  eomputing  the  exaet  solution  of  the  eontinuous  governing  equations.  In  praetiee,  the  TE  is 
approximated  and  the  eomputed  solution  should  eonverge  towards  the  eontinuous  solution  at  a  higher 
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order  rate.  DC  methods  are  extremely  simple  to  implement  as  they  only  require  the  formulation  of  the  TE 
estimate  and  the  ability  to  inelude  a  souree  term  in  the  diserete  solver.  In  addition,  DC  methods  are 
generally  mueh  less  eostly  to  solve  than  the  original  diserete  system  as  they  ean  be  initialized  using  the 
already  available  diserete  solution. 

For  the  results  given  below,  the  quasi- ID  Euler  equations  are  solved  using  a  MUSCL  seheme  (van 
Leer,  1979)  with  k  =  1/3  (eorresponding  to  parabolie  reeonstruetions  on  a  Cartesian  mesh),  the  van 
Albada  flux  limiter  (van  Albada  et  ah,  1982),  and  Roe’s  approximate  Riemann  solver  flux  seheme  (Roe, 
1981).  The  Maeh  number  distributions  and  nozzle  geometry  are  shown  in  Figure  19  for  inflow  stagnation 
eonditions  of  po  =  300  kPa  and  To  =  600  K.  The  nozzle  is  deseribed  by  a  Gaussian  bump  that  avoids 
eurvature  diseontinuities  that  ean  eause  TE  spikes: 


A(x) 


1.0 -0.8  exp 


1  ^ 

1  2  JJ 

-l<v<l.. 


(64) 


Figure  19.  Quasi- ID  nozzle  area  and  Maeh  number  distributions. 

The  funetional  of  interest  for  the  adjoint  error  estimation  is  the  integral  of  pressure  along  the  nozzle 
normalized  by  the  inflow  stagnation  pressure  as  given  by: 


and  has  a  value  of  1.0187219  for  the  isentropie  ease  and  1.2972908  for  the  shoeked  ease  as  ealeulated 
using  a  seven  point  Gauss  quadrature  of  the  analytie  solution  on  a  grid  with  approximately  1  million  eells. 
For  the  diserete  value  of  the  integral,  the  trapezoidal  rule  is  applied  to  the  reeonstrueted  solution. 
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Figure  20a  shows  the  behavior  of  the  base  solution  error  and  remaining  error  after  adjoint  eorreetion 
for  the  funetional  using  the  exaet  TE,  the  weak  estimate  of  the  TE,  the  embedded  grid  approaeh,  and  the 
eorreeted  embedded  grid  approaeh  on  a  series  of  Cartesian  meshes.  As  expeeted,  the  base  error  eonverges 
at  a  seeond  order  rate.  The  standard  embedded  grid  approaeh  only  demonstrates  a  half  order  of  magnitude 


Base  DE 
DC;  exact  TE 


(a) 


(b) 


(c) 


(d) 


Figure  20.  Funetional  base  error  and  remaining  error  for  isentropie  expansion  for  the  adjoint  method 
with  various  TE  estimation  teehniques  (a)  and  for  the  adjoint,  DC,  and  ETEs  in  (b).  Note  the 
differenee  in  seales.  Demonstration  of  sign  ehange  in  TE  eausing  TE  spike  (e)  and  funetional  base 
error  and  remaining  error  for  eropped  domain  whieh  removes  the  TE  spike  (d). 

reduetion  in  the  remaining  error  with  a  seeond  order  eonvergenee  rate,  a  result  seen  elsewhere  in  literature 
(Nemee  and  Aftosmis,  2007).  For  the  exaet  TE  and  weak  form  estimate,  the  remaining  error  eonverges  at 
an  approximately  third  to  fourth  order  rate,  before  beginning  to  level  off  towards  a  first  to  seeond  order 
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rate  at  finer  grid  levels.  The  eorreeted  embedded  grid  method  shows  a  third  order  eonvergenee  rate  and 
surpasses  the  weak  form  TE  estimate  in  a  region  where  the  weak  form  begins  to  demonstrate  first  order 
eonvergenee.  This  lower  order  eonvergenee  is  due  to  a  sign  ehange  in  the  TE  whieh  ereates  a  zeroth  order 
spike  in  the  TE,  as  shown  in  Figure  20e  near  x  =  0.95m.  The  diserete  residual  ineludes  a  limiter  funetion 
whieh  smooths  out  this  oseillation  in  the  embedded  grid  TE  estimate  and  prevents  it  from  dominating  the 
error  estimate.  If  the  domain  is  shortened  to  remove  the  region  whieh  eontains  the  TE  spike,  then  the 
performanee  of  the  weak  form  based  adjoint  error  estimate  is  improved,  as  shown  in  Figure  20d,  but  it 
now  demonstrates  a  sealloped  region  due  to  a  ehange  in  the  sign  of  the  predieted  remaining  error.  It 
should  be  noted  that  the  results  of  Figure  20d  are  obtained  through  new  solutions  of  the  primal  and  dual 
problems  on  the  shortened  domain  instead  of  simply  removing  the  eontributions  of  the  adjoint  weighted 
TE  from  the  eropped  region  of  the  domain. 

Comparisons  between  the  ETEs,  DC,  and  the  adjoint  method  are  made  in  Figure  20b.  Most  notably, 
when  the  exaet  TE  is  used,  the  DC  essentially  returns  the  exaet  funetional  value,  and  the  ETEs  and  adjoint 
method  eonverge  exaetly  the  same.  The  better  performanee  of  the  DC  method  ean  be  attributed  to  the  faet 
that  it  is  solving  the  nonlinear  governing  equations,  while  both  the  adjoint  and  ETEs  are  solving 
linearized  problems.  A  lower  order  approximation  for  the  ETEs  is  ineluded  where  the  linearized  equations 
only  eonsists  of  first  order  Jaeobians.  On  eoarser  meshes,  the  low  order  formulation  of  the  ETEs  will  not 
eonverge  and  perform  about  two  orders  of  magnitude  worse  than  the  fully  linearized  ETEs  on  the  finest 
mesh.  When  an  approximate  TE  estimate  is  used,  in  this  ease  it  is  the  weak  formulation  diseussed  above, 
the  ETEs,  DC,  and  the  adjoint  method  all  perform  similarly  well,  but  the  DC  does  outperform  the  others 
on  the  eoarser  mesh  levels.  Given  the  similarity  in  performanee,  it  appears  that  either  the  DC  method  or 
ETEs  eould  be  used  in  plaee  of  an  adjoint  method  in  order  to  provide  funetional  error  estimates.  Results 
for  the  first  order  ETEs  are  not  ineluded  for  the  approximate  TE  method  as  the  DE  estimates  were  mueh 
worse  than  the  rest. 


X.  Mesh  Adaptation 
Introduction 

Performing  mesh  adaptation  involves  two  distinet  aspeets:  ehoosing  a  method  for  driving  the 
adaptation  and  ehoosing  a  meehanism  to  aetually  perform  the  adaptation.  While  this  report  is  elearly 
foeused  on  the  former,  here  we  briefly  mention  the  meehanism  for  performing  the  adaptation  that  was 
used.  As  diseussed  before,  we  limited  the  seope  of  this  effort  by  demonstrating  our  error  eontrol 
teehniques  on  struetured  meshes.  Sinee  h-adaptation  is  not  possible  on  struetured  meshes,  r-adaptation 
(i.e.,  mesh  movement)  was  used.  Simple  spring  analogy  and  elliptieal  adaptation  methods  (Baker,  1997) 
ean  be  used  when  only  a  single  adaptation  driver  funetion  is  available  (e.g.,  the  total  truneation  error). 
When  more  detailed  information  is  available  regarding  mesh  quality  eontributions,  then  the  spring 
analogy  approaeh  ean  be  supplemented  with  torsional  springs  to  mitigate  skewness  effeets  (e.g., 
Nakahashi  and  Deiwert,  1986).  Furthermore,  variational  methods  allow  for  mesh  adaptation  with  the 
addition  of  mesh  quality  eonstraints  (e.g.,  Lapenta,  2003). 

Our  approaeh  for  driving  the  mesh  adaptation  is  foeused  on  the  residual/truneation  error,  ineluding 
adjoint  weighting  when  speeifie  solution  funetionals  are  of  interest.  We  maintain  that  the 
residual/truneation  error  is  the  only  theoretieally  sound  approaeh  to  driving  mesh  adaptation.  Adaptation 
based  on  solution  features  has  been  shown  to  fail  (Dwight,  2008),  sometimes  “disastrously”  (Ainsworth 
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and  Oden,  2000).  Adaptation  based  on  the  diseretization  error  negleets  the  effeets  of  error  transport  from 
other  regions  of  the  domain  (e.g.,  through  eonveetion  and  diffusion)  and  has  also  been  shown  to  be  a  poor 
driver  for  mesh  adaptation  (e.g.,  Gu  and  Shih,  2001).  Adaptation  based  on  reeovery-based  error  estimates 
(e.g.,  Zienkiewiez  and  Zhu,  1992)  has  seen  some  sueeess  for  linear  elliptie  problems  in  finite  elements, 
but  this  sueeess  is  not  expeeted  to  earry  over  to  hyperbolie  problems  using  finite  differenee  or  finite 
volume  methods. 

When  the  truneation  error  is  available  (e.g.,  from  a  eontinuous  or  diserete  residual),  then  it  ean  be 
used  direetly  to  drive  the  adaptation  proeess.  We  explored  different  methods  for  assembling  the  truneation 
error  for  eaeh  eonservation  equation  (mass,  momentum,  and  energy)  into  a  single  weighting  funetion  to 
drive  adaptation.  Sinee  eaeh  eonservation  equation  has  different  magnitudes  (e.g.,  mass  equation  versus 
total  energy  equation),  we  investigated  methods  for  sealing  the  truneation  error  based  on  appropriate 
referenee  truneation  error  values  from  a  baseline  (i.e.,  unadapted)  mesh.  Onee  this  weight  funetion  is 
assembled,  appropriate  smoothing  or  filtering  methods  are  used  to  provide  for  smooth  adaptation. 

Limiting  of  the  weighting  funetion  is  also  needed  to  prevent  small  values  (whieh  ean  oeeur  when  the 
truneation  error  goes  from  positive  to  negative)  from  affeeting  the  adaptation  proeess.  For  example, 
spring-based  adaptation  when  two  nodes  are  eonneeted  with  a  spring  with  near  zero  spring  eonstant  ean 
result  in  large  loeal  mesh  deformations  and  result  in  mesh  erossover.  When  one  is  instead  interested  in 
targeted  mesh  adaptation  (i.e.,  mesh  adaptation  to  reduee  the  error  in  a  ehosen  solution  funetional),  then 
the  adjoint  solution  provides  the  appropriate  weighting  for  how  the  loeal  eell’s  residual/truneation  error 
will  affeet  the  funetional.  For  example,  instead  of  adapting  the  mesh  purely  on  the  magnitude  of  the 
truneation  error,  one  might  instead  adapt  the  mesh  based  on  the  truneation  error  weighted  by  the  adjoint 
solution. 


When  eondueting  adaptation,  it  is  important  to  monitor  the  eonvergenee  of  the  adaptation  seheme  to 
know  when  to  stop  adaptation.  However,  an  appropriate  adaptation  monitor  must  be  seleeted  that  is 
eonsistent  with  the  adaptation  sehemes  being  used  in  order  to  do  this.  For  this  investigation,  two 
adaptation  monitors  are  seleeted  based  upon  the  eoneept  of  equidistribution.  The  first  adaptation  monitor 
is  a  truneation  error  metrie  given  by  Equation  (66)  whieh  eompares  the  L2  norm  of  truneation  error  on  the 
adapted  grid  to  the  L2  norm  of  truneation  error  on  the  initial  grid.  The  L2  norm  of  truneation  error  in  2D  is 
given  by  Equation  (67)  where  p  is  equal  to  2,  represents  the  truneation  error  in  a  given  eell,  Aq  is  the 
area  of  a  given  eell,  and  A  is  the  total  area  of  the  domain. 
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The  seeond  adaptation  monitor  is  the  weight  funetion  equidistribution  error  given  by  Equation  (68) 
and  is  taken  from  previous  work  by  Choudhary  (2014): 


WEE.j  = 


^-WuAj 


(68) 
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where  a  is  the  equidistributed  value  of  the  quantity  weight  funetion  times  eell  area  aeross  the  entire 
eomputational  domain.  For  all  test  eases  in  this  study,  the  eonvergenee  toleranee  for  these  adaptation 
monitors  was  set  to  0.1.  Beyond  the  toleranee  of  0.1,  mesh  movement  is  minimal  and  the  benefit  of 
further  adaptation  does  not  outweigh  the  added  eomputational  eost. 


Sinee  truneation  error  is  used  as  the  driver  for  adaptation,  it  is  neeessary  to  formulate  a  weight 
funetion  in  terms  of  truneation  error.  To  holistieally  adapt  the  mesh  for  the  entire  system,  the  truneation 
error  from  eaeh  governing  equation  must  be  eombined  in  a  eonsistent  manner.  This  is  aeeomplished  by 
defining  the  loeal  weight  funetion  as  the  average  of  the  absolute  value  of  eaeh  equation’s  truneation  error 
normalized  by  its  maximum  absolute  value  on  the  initial  or  uniform  grid.  This  is  given  by  Equation  (69) 
where  neq  is  the  number  of  governing  equations  in  the  system: 


K\ 

max 

\te'‘ 

(69) 


When  eonstrueting  the  weight  funetion  in  this  fashion,  often  there  is  some  amount  of  high  frequeney 
noise  present.  To  prevent  the  adaptation  sehemes  from  diverging  or  eausing  mesh  erossover,  it  is 
neeessary  to  apply  a  smoothing  algorithm  to  the  weight  funetion  prior  to  its  applieation.  Weight  funetion 
smoothing  is  aeeomplished  by  applying  the  elliptie  smoother  given  in  Equation  (70): 


'^smooth  _  i-i,J  i+lj  hj-l  hj+l  i,J 


12 


(70) 


Although  some  weight  funetion  smoothing  is  required,  it  is  also  important  that  the  weight  funetion  is  not 
oversmoothed.  If  the  weight  funetion  is  over-smoothed,  information  about  the  problem  is  lost,  and  the 
amount  of  diseretization  error  improvement  will  likely  be  diminished. 


Adaptation  Results 

This  report  examines  the  effeetiveness  of  several  2-D  r-adaptation  sehemes  in  redueing  diseretization 
error  in  numerieal  solutions.  The  adaptation  sehemes  used  inelude  an  adaptive  Poisson  grid  generator 
(Anderson,  1990),  a  variational  grid  generator  (Braekbill  and  Saltzman,  1982),  a  eenter  of  mass  based 
seheme  (Lafiin,  1997),  and  a  seheme  based  on  deforming  maps  (Liao  and  Anderson,  1992).  These 
sehemes  are  applied  to  Maeh  1.2  flow  around  a  12  deg.  downward  turn.  Diseretization  error  is  eomputed 
using  the  known  exaet  solution.  Diseretization  error  reduetions  of  up  to  thirteen  times  are  aehieved  on 
adapted  grids  relative  to  the  diseretization  error  present  on  a  uniform  grid  of  the  same  size.  This  error 
reduetion  results  in  a  high  degree  of  effieieney  sinee  the  primal  problem  eonverges  at  a  sublinear  rate  (i.e., 
less  than  first  order)  as  expeeted  (Banks  et  al.,  2008). 

The  test  problem  is  a  2D  supersonie  expansion  fan  generated  by  a  Maeh  1.2  flow  around  a  12° 
downward  turn.  Three  instanees  of  this  problem  are  examined  and  will  be  denoted  as  the  “edge  ease,”  the 
“eomer  ease,”  and  the  “shifted  ease.”  The  eomputational  domain  for  these  eases  ean  be  found  in  Figure 
21 .  The  grid  for  eaeh  ease  setup  is  a  65x65  node  grid  extending  from  0  m  to  1  m  in  both  the  x  and  y 
direetions.  The  inflow  eonditions  for  all  eases  are  statie  pressure  and  temperature  of  100  kPa  and  273  K, 
respeetively. 
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Figure  21.  Supersonic  expansion  fan  domains:  black  square  denotes  the  corner  case,  blue  square 
denotes  the  edge  case,  and  the  green  square  denotes  the  shifted  case. 


Edge  Case 

Truncation  error  and  discretization  error  improvements  for  the  edge  case  of  the  supersonic  expansion 
fan  achieved  by  each  adaptation  scheme  may  be  found  in  Tyson  et  al.  (2015).  In  the  best  case,  truncation 
error  is  decreased  by  a  factor  of  ten  when  using  the  center  of  mass  approach  (Laflin,  1997).  For  the  same 
adaptation  scheme,  discretization  error  is  reduced  by  a  factor  of  six.  Anderson’s  adaptation  scheme 
(Anderson,  1990)  also  performs  well  by  achieving  a  six  time  improvement  in  discretization  error.  The 
Brackbill  and  Saltzman  (1982)  and  Deformation  (Liao  and  Anderson,  1992)  adaptation  schemes  still 
achieve  some  amount  of  discretization  error  reduction,  but  at  best  it  is  half  the  reduction  seen  with  the 
other  two  methods. 

Qualitative  results  of  truncation  error  and  discretization  error  reductions  for  the  edge  case  of  the 
supersonic  expansion  fan  are  presented  in  Figure  22  and  Figure  23,  respectively.  Truncation  error  results 
are  only  presented  for  the  x-momentum  equation  while  discretization  error  results  are  only  presented  for 
pressure.  For  this  case,  the  majority  of  truncation  error  is  located  at  the  beginning  and  end  of  the 
expansion  fan,  as  seen  in  Figure  22e.  Since  truncation  error  is  high  in  these  regions  of  the  domain,  the 
adaptation  schemes  pull  nodes  to  the  beginning  and  end  of  the  expansion  fan  as  illustrated  by  the  adapted 
grids  in  Figure  22a  -  Figure  22d. 
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Figure  22.  Truncation  error  for  the  x-momentum  equation;  edge  case. 
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Figure  23.  Discretization  error  for  the  pressure:  edge  case. 
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Discretization  error  comparisons  for  u-velocity  along  a  line  through  the  domain  at  x  =  0.5  m  is  shown 
in  Figure  24.  The  spikes  in  discretization  error  in  this  figure  at  approximately  y  =  0.05  m  and  y  =  0.95  m 
correspond  to  the  edges  of  the  expansion  fan.  It  can  be  seen  that  the  Anderson  and  center  of  mass 
adaptation  schemes  best  reduce  these  spikes  in  discretization  error.  Figure  25  illustrates  the  adaptive 
convergence  of  the  Anderson  adaptation  scheme  for  this  case.  The  L2  norm  of  truncation  error  and  the 
weight  function  equidistribution  error  are  plotted  against  the  number  of  adaptation  cycles.  Although  it 
may  appear  that  these  monitors  could  converge  further,  truncation  error  norms  do  not  continue  to 
converge  and  no  increased  benefit  is  found  by  adapting  more  than  50  adaptation  cycles  for  this  case. 


Figure  24.  Discretization  error  for  w-velocity  dtx  =  0.5  m:  edge  case. 
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1 


Figure  25.  Adaptive  eonvergenee:  Anderson  method  (edge  ease). 


Corner  Case 

Truneation  and  diseretization  error  improvements  for  the  eorner  ease  of  the  supersonie  expansion  fan 
aehieved  by  eaeh  adaptation  seheme  may  be  found  in  Tyson  et  al.  (2015).  Similar  to  the  edge  ease,  the 
greatest  truneation  error  reduetion  is  aehieved  using  the  eenter  of  mass  approaeh  with  about  an  eight  time 
reduetion.  Diseretization  error  for  this  adaptation  seheme  is  redueed  by  a  faetor  of  approximately  five  for 
eaeh  primitive  variable.  Although  Anderson’s  adaptation  seheme  does  not  perform  as  well  here  as  with 
the  edge  ease,  over  a  three  time  improvement  in  diseretization  error  is  aehieved.  The  Braekbill  and 
Saltzman  and  Deformation  adaptation  sehemes  aehieve  a  peak  reduetion  faetor  in  diseretization  error  of 
1.87  and  1.49,  respeetively. 

Most  of  the  diseretization  error  for  this  ease  is  present  in  the  bottom  left  eorner  of  the  domain  at  the 
start  of  the  expansion  fan.  With  this  in  mind,  diseretization  error  is  only  plotted  for  eaeh  adaptation 
seheme  in  the  bottom  left  eorner  of  the  domain  to  better  illustrate  the  results  for  eaeh  adaptation  method. 
Diseretization  error  in  density  is  given  in  Figure  26.  It  ean  be  seen  in  Figure  26a  and  Figure  26e  how  well 
the  Anderson  and  eenter  of  mass  sehemes  reduee  diseretization  error  in  this  region  relative  to  the 
Braekbill  and  Saltzman  and  Deformation  adaptation  sehemes.  Figure  26b  and  Figure  26d,  respeetively. 
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c)  Center  of  Mass 
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Figure  26.  Discretization  error  for  the  density:  corner  case. 
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Figure  27.  Discretization  error  for  pressure  at  x  =  0.5  m:  corner  case. 


Figure  10.  Discretization  Error  in  Pressure  at  x  =  0.5m 


Figure  28.  Adaptive  convergence:  center  of  mass  method  (corner  case). 
Figure  11.  Adaptive  Convergence:  Center  of  Mass 
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Discretization  error  comparisons  for  pressure  along  a  line  through  the  domain  at  x  =  0.5m  may  be 
found  in  Figure  27.  Similar  to  the  edge  case,  the  spikes  in  discretization  error  in  this  figure  at 
approximately  y  =  0.25m  and  y  =  0.75m  correspond  to  the  edges  of  the  expansion  fan.  Again,  the 
Anderson  and  center  of  mass  adaptation  schemes  are  able  to  reduce  these  spikes  in  discretization  error  the 
most.  Figure  28  illustrates  the  adaptive  convergence  of  the  center  of  mass  adaptation  scheme  for  this  case. 
The  L2  norm  of  truncation  error  and  the  weight  function  equidistribution  error  are  plotted  against  the 
number  of  adaptation  cycles.  Here  adaptation  is  performed  for  150  adaptation  cycles.  Past  50  adaptation 
cycles,  the  adaptation  is  converged  and  no  added  benefit  is  achieved  by  further  adaptation. 

Shifted  Case 

The  discretization  error  in  pressure  for  the  shifted  case  is  shown  in  Figure  29  for  the  uniform  mesh 
and  the  center  of  mass  adaptation  method.  The  center  of  mass  approach  achieved  a  factor  of  13.5 
reduction  in  the  L2  norm  of  the  discretization  error  in  the  pressure.  These  higher  improvement  factors 
relative  to  the  corner  case  are  due  to  the  singularity  being  located  on  only  one  boundary  instead  of  the 
intersection  of  two  boundaries.  The  discretization  error  in  the  pressure  extracted  from  a  line  at  x  =  0.5  m 
is  shown  in  Figure  30.  Large  reductions  in  the  error  are  found  at  the  beginning  and  ending  of  the 
expansion  fan  relative  to  the  uniform  mesh. 


X 


Figure  29.  Discretization  error  in  pressure  for  the  shifted  case:  uniform  mesh  (left)  and  center  of  mass 

adapted  mesh  (right). 
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Figure  30.  Discretization  error  for  pressure  at  jc  =  0.5  m:  shifted  case. 


Grid  Study 

A  grid  study  is  conducted  for  the  corner  case  of  the  2D  expansion  fan.  The  L2  norm  of  discretization 
error  in  pressure  is  compared  on  uniform  and  adapted  grids  for  five  uniformly  refined  grid  levels  ranging 
from  9x9  to  129x129.  Each  adapted  grid  in  this  study  is  generated  independently  starting  from  a  uniform 
grid.  Figure  31  illustrates  the  savings  which  can  be  achieved  when  using  r-adaptation  as  a  means  of 
reducing  discretization  error.  For  example,  to  achieve  the  same  level  of  discretization  error  as  the  coarsest 
adapted  grid,  a  uniform  grid  which  is  three  grid  levels  finer  would  be  required.  This  uniform  grid  would 
have  64  times  the  number  of  cells  of  the  coarsest  adapted  grid  for  this  2D  problem.  It  would  be  512  times 
as  many  cells  in  3D.  This  equates  to  huge  computational  savings  especially  when  a  lower  level  of 
discretization  error  is  desired.  Similar  trends  can  also  be  seen  with  the  other  primitive  variables  for  this 
case.  Also,  the  order  of  convergence  between  the  uniform  and  adapted  grids  is  quite  similar.  This  is 
surprising  considering  the  adaptation  schemes  used  in  this  study  provide  no  guarantee  that  the  order  of 
convergence  on  adapted  grids  will  be  the  same  as  consistently  refined  unadapted  grids.  The  order  of 
convergence  using  discretization  error  in  pressure  for  the  three  finest  grid  levels  is  computed  to  be 
p  =  0.68 .  This  convergence  rate  is  consistent  with  what  is  expected  for  problems  containing  linear 
discontinuities  such  as  this  expansion  fan  case.  Banks  et  al.  (2008)  show  that  for  problems  with  linear 
discontinuities  the  observed  order  of  accuracy  will  reduce  to: 


P  = 


Pf 

Pf+i 


(71) 


where  pf  is  the  formal  order  of  accuracy  of  the  discretization.  With  a  formal  order  of  accuracy  of  pf=  2, 
the  expected  order  of  accuracy  for  this  case  will  be  p  =  0.66 . 
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Figure  3 1 .  Grid  study  showing  the  behavior  of  the  L2  norm  of  diseretization  error  in  pressure  for 
uniform  and  eenter  of  mass  adapted  meshes:  eorner  ease. 
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